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Preface 


The  purpose  of  this  research  has  been  to  revisit  the  atmospheric  nuclear  test  data  of  the 
1950’s  and  1960’s  in  order  to  gain  information  germane  to  nuclear  winter.  I  have  concen¬ 
trated  on  resolving  the  size  distribution  of  the  debris  ensemble  because  particle  size  strongly 
governs  the  magnitude  and  duration  of  sunlight  attenuation.  By  modeling  the  removal  of 
nuclear  debris  from  the  atmosphere,  it  has  been  possible  to  bound  the  admissible  distribution 
of  lofted  particulates.  In  addition,  I  have  shown  how  gravitational  cloud  stratification 
creates  particle  populations  approximated  by  a  power  law  distribution.  The  results  should 
be  useful  for  fallout  and  dust  effects  modeling  as  well  as  nuclear  winter  studies.  I  was  frus¬ 
trated  by  time  limitations—  there  is  much  more  to  be  learned  from  the  '400  events  which 
have  occurred  in  our  atmosphere. 
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Marcel  Nathans  (LLNL),  Dr.  Ernest  Bauer  (IDA),  and  Dr.  Alan  Mason  (LASL)  provided  the 
foundation  for  much  of  what  I  have  done.  1  am  most  grateful  for  their  advice  and  counsel 
during  several  meetings  and  telephone  conversations.  Capt  Norm  Davis  (USMC)  provided 
valuable  assistance  during  the  course  of  his  masters’  research.  In  addition,  without  the 
material  provided  by  Mr.  A1  West  of  the  DNA  Technical  Library  and  Mssrs.  Bill  Alfont  and 
Dick  Rowland  of  KAMAN  TEMPO  my  research  would  not  have  been  possible.  The  interest 
and  assistance  of  LTC  Ron  Tuttle  of  AFIT  and  Dr.  David  Auton  of  the  Defense  Nuclear 
Agency  was  most  appreciated.  Special  thanks  are  due  to  Dr.  Marvin  Atkins,  Dr.  Gordon 
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my  degree. 
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ABSTRACT 


Atmospheric  test  fallout  data  have  been  used  to  determine  admissible  dust 
particle  size  distributions  for  nuclear  winter  studies.  The  research  was  originally 
motivated  by  extreme  differences  noted  in  the  magnitude  and  longevity  of  dust 
effects  predicted  by  particle  size  distributions*  routinely  used  in  fallout  predictions 
versus  those  used  for  nuclear  winter  studies.  Three  different  sets  of  historical  data 
have  been  analyzed: 

1)  Stratospheric  burden  of  Strontium-90  and  Tungsten-185,  1954-1967 
(97  contributing  events)  ' 

2)  Continental  U.S  Strontium-90  fallout  through  1958 
(75  contributing  events)  ^  ar,-,d^ 

h)  Local  Fallout  from  selected  Nevada  tests  (16  events)  4 


The  contribution  of  dust  to  possible  long  term  climate  effects  following  a 
nuclear  exchange  depends  strongly  on  the  particle  size  distribution.*  The  distribu¬ 
tion  affects  both  the  atmospheric  residence  time  and  optical  depth.  One  dimen¬ 
sional  models  of  stratospheric/tropospheric  fallout  removal  were  developed  and 
used  to  identify  optimum  particle  distributions.  Results  indicate  that  particle  dis¬ 
tributions  which  properly  predict  bulk  stratospheric  activity  transfer  tend  to  be 
somewhat  smaller  than  number  size  distributions  used  in  initial  nuclear  winter 
studies.  In  addition,  both  90 Sr  and  185  W  fallout  behavior  is  better  predicted  by 
the  lognohnal  distribution  function  than  the  prevalent  power  law  hybrid  function. 
It  is  shown  thjiKthe  power  law  behavior  of  particle  samples  may  well  be  an  aber¬ 


ration  of  gravitationa^cloud  stratification.  Results  support  the  possible  existence 
of  two  independent  particle  size  distributions  in  clouds  generated  by  surface  or 
near  surface  bursts.  One  distribution  governs  late  time  stratospheric  fallout,  the 

v... 

other  governs  early  time  fallout.  A'bimodal  lognormal  distribution  is  proposed  to 
describe  the  cloud  particle  population.  'T^e  distribution  predicts  higher  initial 
sunlight  attenuation  and  lower  late  time  attebvjation  than  the  power  law  hybrid 
function  used  in  initial  nuclear  winter  studies.  v 
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IMPLICATIONS  OF  ATMOSPHERIC  TEST  FALLOUT 
DATA  FOR  NUCLEAR  WINTER 


I.  Introduction. 


A.  Background. 

The  possibility  of  serious  climatic  consequences  of  a  nuclear  war  has  been 
recently  investigated  by  several  groups  (1,2, 3, 4, 5).  Perhaps  the  best  known  investiga¬ 
tion,  and  one  with  the  most  cataclysmic  predictions,  was  by  R.  P.  Turco,  O.  B. 
Toon,  TP.  Ackerman,  JP.  Pollack,  and  C.  Sagan.  The  "TTAPS"  group  concluded 
that  attenuation  of  incoming  solar  radiation  leading  to  subfreezing  land  temperatures 
oould  occur  over  the  northern  hemisphere  due  to  vast  amounts  of  dust  and  smoke 
aerosols  injected  into  the  atmosphere  during  a  major  nuclear  exchange.  Their  initial 
calculations  (1,  see  figure  1  inset)  were  based  upon  a  9800  megaton  baseline  exchange 
between  the  superpowers  but  indicated  that  the  yield  threshold  for  major  optical  (and 
therefore  climatic)  effects  may  be  as  low  as  100  megatons.  The  outcome  the  TTAPS 
baseline  scenario  is  a  peak  35°  C  temperature  drop  of  roughly  6  month  duration 
(figure  2).  A  subsequent,  independent  evaluation  of  these  findings  by  the  National 
Research  Council  (3)  questioned  the  plausibility  of  certain  of  the  TTAPS  assumptions 
(e.g.  treatment  of  smoke  injection),  but  confirmed  the  possibility  of  large  and  pro¬ 
longed  temperature  drops  across  the  northern  hemisphere  leading  to  wbat  the  TTAPS 
group  dubbed  "nuclear  winter."  Ramaswamy  and  Kiehl  (abbreviated  R-K),  in  a 
recent  parametric  evaluation  of  the  optical  properties  of  dust  and  smoke,  also  calcu¬ 
lated  significant  temperature  drops  (5).  In  these  initial  studies,  smoke  appears  to 
dominate  dust  because  of  its  higher  absorption  properties,  and  the  large  fuel  deusities 
assumed  for  targeted  urban  areas.  Thus  most  of  the  research  to  date  has  been 
devoted  to  verifying  smoke  effects.  This  lopsided  attention  to  smoke  has  also  been 
warranted  by  the  large-  number  of  variables  in  the  fire  problem  which  has  left  much 
more  room  for  interpretation  and  hence  much  larger  uncertainties.  Differing 


1 


TYPE  OP  BURST 


NUMBER 

Of 

BURSTS 


YIELD  PER 
WARHEAD 
(megatons) 


YIELD 

EXPENDED 

(megatons) 


Land  surface 
Land  surface 
Land  near' surface 
Low  alrburat 
Exo- atmospheric 
Low  airburst 
Low  airburst 
Low  airburst 
Water  surfaee 
Low  alrburat 
Water  surface 


10 

S 

1 

1 

1 

O.B 

0.3 

0.2 

0.2 

0.1 

0.1 


110 

400 

3000 

1000 

100 

500 

3000 

3000 

1000 

3000 

1000 

10100 


1100 

2200 

3000 

1000 

100 

250 

900 

000 

200 

300 

100 

9000 


Figure  1.  TTAPS  Nuclear  Exchange 


assumptions  concerning  whether  fires  start  at  all,  whether  they  spread  significantly, 
whether  firestorms  occur  and,  if  so,  whether  the  smoke  reaches  the  stratosphere  have 
led  to  extreme  variations  in  the  severity  of  predicted  outcomes. 

While  it  may  be  true  that  the  amplitude  of  dust  induced  climate  effects  is  less 
severe,  the  longevity  of  stratospheric  dust  aerosols  (as  compared  to  longevity  of 
smoke  which  is  confined  predominantly  to  the  troposphere  and  is  therefore  subject  to 
rapid  removal  by  washout)  prolongs  the  dust  induced  temperature  drop  for  periods 
exceeding  one  year.  Indeed,  dust  alone  is  sufficient  to  produce  significant  temperature 
changes  according  to  the  results  of  a  TTAPS  scenario  in  which  no  smoke  was  injected 
(figure  2).  The  relative  significance  of  dust  has  been  enhanced  by  declining  estimates 
of  urban  fuel  loading  and  the  addition  of  scavenging  to  smoke  removal  models  (6,7). 

Although  not  as  formidable  as  those  for  smoke,  large  uncertainties  do  exist  in 
the  parameters  used  for  the  modeling  of  dust  effects.  One  of  the  largest  of  these 
uncertainties  is  the  particle  size  distribution.  The  optical  thickness  of  a  given  mass  of 
dust  is  extremely  sensitive  to  the  size  distribution  of  the  dust.  This  is  illustrated  in 
Chapter  II  where  the  optical  thickness  is  computed  for  three  different  size  distribu¬ 
tions  as  proposed  by  Nathans  (and  used  by  TTAPS),  Ramaswamy  (5),  and  Norment 
(8).  The  large  differences  in  the  magnitude  and  duration  of  the  sunlight  occlusion  for 
the  three  selected  size  distributions  and  the  obvious  extreme  sensitivity  of  optical 
thickness  behavior  to  the  dispersion  of  the  size  distribution  are  worrisome  in  that  the 
relative  effect  on  the  climate  ranges  from  nil  (Norment)  to  potentially  severe 
(Nathans,  Ramaswamy). 

B.  Research  Objective. 

The  objective  of  the  present  effort  has  been  to  reduce  the  Uncertainties  in  the 
admissible  weapon  debris  aerosol  size  distributions.  While  large  uncertainties  still 
remain  with  respect  to  the  mass  lofted  vs.  yield  and  burst  height,  the  uncertainties  in 
size  distribution  are  more  significant  because  they  affect  both  the  amplitude  and 
duration  of  sunlight  attenuation.  Variance  in  the  lofted  mass  affects  only  the  ampli¬ 
tude  of  the  attenuation  and  probably  does  not  exceed  an  order  of  magnitude  for  sur¬ 
face  bursts  (1,3).  This  translates  into  a  factor  of  ten  uncertainty  in  the  optical  thick¬ 
ness.  It  is  evident  from  figure  3  that  larger  uncertainties  in  the  optical  thickness  can 
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Figure  2.  Dust  vs  Smoke  Effects 


be  attributed  to  variance  in  size  distribution. 


C.  Approach. 

This  effort  has  used  activity  tracing  to  bound  admissible  size  distributions. 
Activity  tracing  involved  starting  with  historical  data  recording  nuclear  weapon  fal¬ 
lout  in  the  atmosphere  and  on  the  ground  and  then  back-calculating  the  size 
distribution(s)  which  could  have  caused  the  measured  activity  transfer.  This  method 
of  determining  admissible  size  distributions  has  not  been  previously  used. 

The  approach  proceeded  as  follows.  First,  fallout  removal/deposition  mechan¬ 
isms  were  modeled.  The  models  were  exercised  against  available  data.  Size  distribu¬ 
tion  parameters  were  adjusted  to  achieve  favorable  comparisons  between  the  models 
and  the  data.  Care  was  taken  to  identify  and  resolve  differences  between  the  size  dis¬ 
tributions  derived  from  activity  tracing,  and  those  previously  reported  (which  were 
based  upon  microscopic  techniques).  A  variety  of  historical  data  was  used.  Chapter 
IV  contains  a  detailed  discussion  of  the  approach,  how  it  differed  from  previous  n(r) 
investigations,  and  how  it  was  affected  by  the  types  of  data  available. 

D.  Other  Methods. 

Other  methods  have  been  used  and  are  being  developed  to  characterize  size  dis¬ 
tributions.  In  the  past,  particle  size  distributions  have  been  measured  primarily  by 
microscopic  examination  of  fallout  samples.  Another  method  presently  under 
development  uses  laser  interferometry.  Both  methods  are  "direct"  in  that  they  meas¬ 
ure  the  size  of  individual  particles  in  isolated  samples.  By  contrast,  the  activity  trac¬ 
ing  method  is  indirect  in  that  size  information  is  inferred  from  bulk  activity  transfer 
data. 


D.l.  Microscopic  Techniques. 

Particle  size  distributions  have  been  primarily  determined  from  the  microscopic 
study  of  individual  samples  from  clouds  (particles  collected  on  aircraft  filters)  or  the 
ground  (gummed  paper,  granular  collectors).  A  discussion  of  aircraft  sampling 
methods  is  contained  in  documentation  of  the  DNA  High  Altitude  Sampling  Program 
(HASP)  (9,10,11,12).  A  discussion  of  ground  fallout  collection  methods  is  contained 


in  test  director  reports  for  many  of  the  atmospheric  test  series  (examples  13,14,15). 


D.l.a.  Nathans’  Cloud  Sample  Analysis. 

Nuclear  winter  studies  have  based  their  particle  size  distributions  on  Nathans’ 
laboratory  analysis  of  cloud  samples  from  four  ground  bursts:  Johnie  Boy,  Castle 
Bravo,  Koon,  and  Zuni  (16).  In  his  experimental  procedure,  sections  of  the  sampling 
filter  paper  were  "ashed"  leaving  behind  the  debris  particles  plus  some  limited  filter 
residue.  Low  temperature  ashing  was  used  when  possible  to  avoid  losing  volatile 
fission  products.  Because  small  particles  were  masked  by  larger  ones  (diameters  range 
over  several  orders  of  magnitude),  samples  were  put  into  solution  and  separated  via 
centrifuge  into  (typically  ten)  size  fractions.  About  100  particles  from  each  size  frac¬ 
tion  were  measured  so  that  the  size  distribution  of  each  sample  was  based  on  '1000 
particles.  By  adding  the  size  distributions  of  the  fractions  weighted  by  number,  com¬ 
posite  size  distributions  were  plotted.  Apparently,  histograms  for  these  surface  bursts 
were  based  upon  irregularly  shaped  particles  only.  Spherical  particles  were  present  in 
significant  numbers  only  in  the  three  smallest  size  fractions  of  Johnie  Boy  and  had  a 
size  distribution  which  differed  from  the  irregulars.  Nathans’  results  for  samples  from 
the  four  surface  bursts  appear  in  figures  4-6.  His  results  include  a  slight  correction 
for  sedimentation  (further  discussion  in  Section  V.E).  The  TTAPS  and  NRC  particle 
size  distributions  were  based  on  these  results. 

Nathans  has  also  done  a  similarly  meticulous  evaluation  of  clouds  from  eleven 
air  bursts  (17,18)  and  finds  in  all  cases  a  unimodal  lognormal  size  distribution  with  a 
slope  similar  to  that  of  the  submicron  surface  burst  cloud  population,  but  having  a 
smaller  median  radius  on  average  (.07/i  air  burst  vs  .25/lx  surface  burst).  Nathans 
claims  to  have  been  able  to  detect  particles  down  to  and  slightly  below  .02 in  diame¬ 
ter.  Below  0.1 /i  he  used  activity  to  isolate  the  particles.  For  extremely  small  parti¬ 
cles  he  developed  a  method  to  agglomerate  several  smaller  particles  in  order  to  detect 
their  presence  (19).  Nathans  noted  a  lower  cutoff  diameter  for  debris  particles  in  the 
vicinity  of  .02 (J.  (18).  Nathans’  analysis  included  the  correlation  of  specific  activity 
(equivalent  fissions  per  gram)  with  particle  radius.  This  correlation  was  an  essential 
input  to  the  activity  tracing  method  pursued  in  the  present  effort. 
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D.l.b.  Ground  Sample  Analysis. 

Unfortunately,  the  documentation  of  n(r)  analysis  for  ground  samples  is  very 
sparse.  Most  analysis  was  concerned  primarily  with  the  radiochemical  morphology  of 
the  particles  and  to  a  lesser  extent  with  the  size  distribution.  In  addition,  reported 
ground  sample  particle  size  distributions  exhibit  much  more  variability  than  do  cloud 
samples.  Norment’s  DELFIC  nominal  distribution,  previously  discussed  (8),  is  a  uni- 
modal  lognormal  function  with  median  radius,  rm  —  .204 ft  and  a  logarithmic  slope 
=  4.  Freiling  (20)  reports  lognormal  distributions  with  mass  median  radii  of  the 
order  of  100f(  and  logarithmic  slopes  of  1.68  to  1.98.  Tompkins  (21)  reports  a  mass 
median  radius  of  250ft  and  a  logarithmic  slope  of  1.9  for  a  Small  Boy  sample  5.8  km 
from  ground  zero.  At  greater  distances,  his  samples  are  not  unimodal.  The  WSEG 
fallout  model  (22),  based  on  ground  deposition,  uses  an  activity  median  of  60ft  and 
slope  of  2.  The  present  research  has  developed  evidence  that  the  extreme  variability 
in  ground  samples  may  well  be  due  to  cloud  stratification  caused  by  gravity  sorting 
(see  section  V.E). 

D. 2.  Laser  Interferometry. 

Optical  methods  involving  laser  interferometry  are  presently  being  perfected  for 
analysis  of  nuclear  cloud  samples  (23).  In  particular,  Los  Alamos  is  developing  a 
method  using  the  doppler  shift  of  laser  light  to  correlate  the  Brownian  motion  of  par¬ 
ticulates  in  solution,  thereby  determining  particle  size  from  diffusion  velocity  (24). 
While  laser  techniques  have  been  successfully  employed  for  in  situ  measurement  of 
particle  sizes  during  recent  high  explosive  events,  systematic  studies  of  old  nuclear 
cloud  samples  by  such  methods  are  not  yet  available.  These  techniques  look  very 
promising. 

E.  Some  Problems  Associated  with  Other  Methods. 

The  advantages  of  the  direct  techniques  are  obvious.  The  analyst  is  measuring 
and  counting  observable  particles  from  real  samples.  The  analyst  knows  when  and 
where  the  sample  originated  as  well  as  which  event  (and  relevant  detonation  parame¬ 
ters)  was  the  primary  contributor  to  observables.  However,  there  are  drawbacks. 
The  size  separation  techniques  may  well  affect  particle  size  (some  dissociation  of 
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agglomerated  particles  may  occur  in  the  solution).  There  are  questions  about  the 
degree  to  which  the  filter  ash  may  mask  the  properties  of  the  debris.  In  addition, 
there  appear  to  be  large  variabilities  between  samples  taken  at  different  times  and 
locations  (even  for  cloud  samples).  The  question  then  remains,  "which  sample  is  to  be 
trusted  as  the  most  truly  representative  of  the  total  cloud?"  The  analysis  procedure 
is  extremely  tedious  and  very  few  events  have  been  systematically  characterized  (five 
surface  burst  clouds,  eleven  air  burst  clouds).  Surface  burst  data  is  of  most  impor- 
tance  for  optical  effects  since  air  bursts  loft  little  mass.  Despite  sample  analysis  for 
many  events,  ground  fallout  histograms  are  so  highly  variable  it  is  difficult  to  con¬ 
struct  a  composite  size  distribution  for  any  one  event  based  upon  microscopic  tech¬ 
niques. 

F.  Purpose. 

As  stated  above  the  purpose  of  this  research  was  not  to  rework  the  direct 
methods  described  above,  but  rather  to  determine  size  distributions  using  a  different 
approach.  This  approach,  activity  tracing,  involved  back-calculating  the  size  distri¬ 
bution  from  recorded  fallout  data.  The  data  and  methods  used  to  trace  activity  are 
presented  in  Chapters  HI  and  IV  respectively. 
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II.  Sensitivity  of  Optical  Thickness  to  Dust  Size 


The  optical  properties  of  any  aerosol  are  extremely  sensitive  to  the  population  vs 
radius.  Assuming  a  uni  modal  lognormal  size  distribution  for  the  aerosol,  it  is  possi¬ 
ble  to  derive  a  relatively  simple  expression  for  the  optical  thickness  of  dust  to  illus¬ 
trate  this  size  sensitivity. 


A.  Optical  Attenuation. 

The  attenuation  of  sunlight  results  from  scattering  and  absorption  by  particu¬ 
lates  lofted  into  the  atmosphere.  If  a  cloud  of  particles  is  separated  into  vertical  lam¬ 
ina,  then  the  differential  change  (due  to  absorption  and  scattering)  in  the  number  of 
photons  traversing  each  lamina  is  related  to  the  aggregate  cross  section  of  the  parti¬ 
cles  in  the  laminar  volume: 

dtp  __  aggregate  particle  arra  nAcds  <fi> 

<p  cloud  cross  section  Ac  ' 


Ae  denotes  cloud  cross  section,  which  is  roughly  equal  to  the  area  of  the  northern 
hemisphere  (1.14  X  10M  m2)  in  a  nuclear  winter  scenario;  n  is  the  particle  number 
density;  and  <fi>  represents  the  average  optical  cross  section  of  the  particles. 
Integration  of  equation  (1)  yields  an  exponential  attenuation: 


R 


e  A'  -  c-°T 


where  OT  is  the  "optical  thickness"  and  Nt  is  the  total  number  of  particles  in  the 
cloud.  Qe  is  the  Mie  extinction  efficiency  which  is  a  function  of  refractive  index  and 
the  ratio  of  diameter  to  wavelength  (25).  The  extinction  efficiency  is  largest  for  sub¬ 
micron  particles  and  averages  about  2.5  at  wavelengths  of  0.5  /xm  for  particle  radii 
between  0.1  and  10  /i  (figure  7).  It  is  evident  from  this  expression  that  the  optical 
thickness  is  directly  proportional  to  the  aggregate  surface  area  of  the  particles  in  the 
cloud.  Thus,  the  rate  of  decay  of  OT  with  time  is  directly  proportional  to  the 
decrease  in  the  second  moment  of  particle  number  distribution  as  particles  are 


removed  from  the  atmosphere  by  gravity  and  diffusion,  i.e., 

00 

OT(t)  a  Jirr*n(r  ,t)dr 
o 

where 


Nt  is  simply: 


J  n(rfi)dr  —  1 
o 


j*<  r3>p. 


where  Mt  is  the  total  mass  lofted,  pt  is  the  debris  density,  and 


<r3>  -fr3n(r,0)dr  (6) 

o 

Assuming  a  lognormal  distribution  (a  discussion  of  the  mathematical  properties  of 
lognormal  distributions  is  provided  in  Appendix  B): 


ri(r,0) 


1  (/n(r)-/n(rm)  f 

^XP  l  2  1  0  J. 


Similarly, 


<r3>  -  exp  |3/n(rm )  4-  —0* 


<r2>  -  exp  [2/n(rm)  +  2/S2] 


and  the  initial  optical  thickness  can  be  expressed  in  terms  of  total  mass  lofted 
(roughly  1012  kg  for  the  5000  megaton  TTAPS  scenario),  effective  particle  radius, 

L- 

cloud  cross  section,  and  particle  density  (2600  — rr  according  to  Nathans,  ref.  16): 


OTtm, 


ZMTQt 

4r'P»Ac 
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where  re  —  rm  exp 


Note  that  the  optical  thickness  of  the  initial  stabilized  cloud  is  dependent  on 
both  the  median  radius  and  dispersion  of  the  size  distribution.  The  ratio  of  the  ini¬ 
tial  optical  thickness  predicted  by  any  two  lognormal  distributions  is  given  by: 


OTx 

ot7 


rM« 

- exp 

fm. 


(11) 


A  mode  radius  of  .204/i  and  a  slope  of  4  is  used  as  the  nominal  parameters  in  the 
Department  of  Defense  Land  Fallout  Interpretive  Code  (DELFIC,  ref.  8). 
Ramaswamy  and  Kiehl  (5)  use  a  pure  lognormal  number  distribution  with  a  median 
radius  of  .16/i  and  a  slope  of  2.  A  factor  of  47  difference  in  optical  thickness  results 
which  translates  to  a  "2  order  of  magnitude  difference  in  predicted  sunlight  attenua¬ 
tion  (using  a  radiative  transfer  algorithm  in  which  full  account  is  taken  of  multiply 
scattered  radiation,  refs.  26,27;  figure  9).  Note  the  exponential  effect  that  the 
assumed  slope  has  on  the  sunlight  attenuation.  As  will  become  evident,  the  slope  also 
has  a  dramatic  effect  on  the  rate  of  dust  removal. 


B.  Subznicron  Particle  Criticality. 


B.l.  Optical  Efficiency. 


As  is  evident  from  figure  7,  the  extinction  efficiency  is  largest  below  1/i  at  the 
peak  solar  wavelength.  Ramaswamy  (5)  illustrates  the  importance  of  submicron  par¬ 
ticulates  to  optical  attenuation  by  computing  the  optica)  thickness  as  the  product  of 
a  mass  attenuation  coefficient,  rpe(vinits  m2/g),  and  a  columnar  particle  area  density: 


OT  -  /0, 


fm(r,z)  dz 


]dr 


(12) 


where  r  is  particle  radius  and  m  is  the  mass  density  of  the  aerosol.  He  then  plots  the 
ipe  as  a  function  of  surface  mode  radius  for  assumed  particle  size  distributions  of 
slope  2.  His  results  are  graphed  in  figure  8  and  vividly  illustrate  the  contribution  of 
the  submicron  particle  population  to  the  optical  thickness.  An  attenuation  maximum 
occurs  at  a  surface  mode  radius  of  0.2/i.  Figure  8  also  illustrates  the  very  small 
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Figure  7.  Extinction  Efficiency  of  Dust  Aerosols  at  Wavelengths 
of  0.5  and  11  pm,  as  a  Function  of  Size 

Source  Journal  of  Qeophyilcal  Research.  Vol.  90,  No.  03  06698.  Copyright  American  Geophysical  Union 


Figure  8.  Dependence  of  the  Specific  Extinction  and  Absorption  of 
(a)  Smoke  Aerosols  ?nd  (b)  Dust  Aerosols  at  a 
Wavelength  of  0.5  pm  on  the  Mode  Radius 


Sourca,  Journal  of  Geophysical  Raaearch,  Vol,  90,  No.  03,  pSSOO,  Copyright  American  Geophysical  Union 
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Figure  9.  Fraction  of  Incident  solar  radiation  reaching  the  surface  as 
a  function  of  extinction  optical  depth  for  smoke  and  dust 
particulates. 
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contribution  of  the  absorption  component,  ,  to  the  attenuation  coefficient  of  dust. 
Dust  is  not  highly  absorbing,  and  the  optical  attenuation  in  dust  clouds  is  primarily 
due  to  scattering.  Thus  for  a  given  optical  thickness,  optical  attenuation  in  dust  is 
much  lower  than  for  smoke  which  it  highly  absorbing  (figure  9). 

B. 2.  Longevity. 

Because  of  the  extremely  low  sedimentation  velocities  of  particles  smaller  than 
lfi,  the  fraction  of  the  particle  ensemble  in  the  submicron  range  governs  the  duration 
of  optical  attenuation  over  periods  necessary  to  induce  climate  changes.  For  instance, 
a  spherical  10m  particle  falls  from  20  km  to  the  tropopause  (average  altitude  12  km) 
in  approximately  3  days  whereas  a  1m  particle  takes  190  days.  A  graph  of  fall  time 
vs  radius  and  altitude  for  spherical  particles  is  included  as  figure  10.  The  graph  is 
based  upon  Stokes-Davies-Macdonald  fall  mechanics  model  as  described  in  Appendix 
A.  It  is  clear  from  this  graph  that  the  submicron  particle  fraction,  in  addition  to 
dominating  the  optical  attenuation,  also  governs  the  duration  of  sunlight  attenua¬ 
tion. 

C.  Ground  vs  Cloud  Sample  Dichotomy. 

Depending  upon  the  origin  of  particle  samples,  size  distributions  determined 
from  the  analysis  of  atmospheric  test  dust  exhibit  markedly  different  submicron  frac¬ 
tions.  Most  available  samples  fall  into  two  major  categories:  (1)  local  fallout  (down 
within  roughly  24  hours  of  event)  ground  samples,  and  (2)  early  time  (within  several 
hours)  cloud  samples. 

The  tendency  has  been  to  use  cloud  samples  as  representative  of  initially  lofted 
particle  distributions.  Intuitively,  one  would  expect  that  because  of  sedimentation, 
ground  samples  would  be  biased  toward  the  larger  end  of  the  size  spectrum  and  cloud 
samples  would  be  oppositely  biased.  This  seems  to  be  true.  However,  one  should 
also  expect  that  since  the  origin  of  both  cloud  and  ground  samples  for  any  event  is 
ultimately  the  same  explosion,  there  should  be  similarities  in  particle  characteristics. 

A  point  in  favor  of  using  ground  sample  distributions  is  that  near  ground  zero, 
fallout  appears  to  result  from  the  vertical  component  of  toroidal  circulation  within 
the  cloud  (as  well  as  sedimentation).  Ground  samples  near  ground  zero  should 
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Figure  10.  Fall  Tima  varsua  Size  and  Altitude 
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therefore  represent  a  homogeneous  mixture  of  lofted  particles,  both  large  and  small 
(28:51).  Indeed,  by  ascribing  ground  sample  size  distributions  to  clouds  and  then 
modeling  deposition  by  sedimentation,  reasonable  predictions  of  ground  fallout  con¬ 
tours  are  obtained  (29,30,31,32).  Ground  samples  cannot  be  easily  dismissed  as  non¬ 
representative  of  the  particle  population.  In  addition,  since  many  cloud  samples  (16) 
were  taken  at  altitudes  considerably  below  the  cloud  vertical  center  (since  for  high 
yield  events,  the  cloud  eluded  the  altitude  capabilities  of  sampling  aircraft),  these 
samples  may  also  be  nonrepresentative  of  the  aggregate  cloud  population. 

There  is  some  additional  cause  for  questioning  the  use  of  distributions  derived 
from  available  analysis  of  surface  burst  cloud  samples.  The  TTAPS  and  NRC  stu¬ 
dies  used  particle  size  distributions  based  upon  surface  burst  cloud  sample  analysis  of 
Marcel  Nathans  (16).  Nathans  measured  and  counted  particles  from  four  surface 
bursts:  Johnie  Boy,  Castle  Bravo,  Koon,  and  Zuni.  Of  these,  the  high  yield  shots 

(of  most  interest  for  nuclear  winter  predictions)  were  detonated  over  coral.  Johnie 
Boy,  the  only  burst  over  silicate  soil,  was  low  yield  (.5  kiloton)  and  slightly  buried. 
Since  the  nuclear  winter  phenomenon  is  based  on  mass  lofted  by  megaton  surface 
bursts  over  silicate  soil  there  is  good  reason  to  question  the  applicability  of  the  meas¬ 
ured  size  distributions.  Admittedly,  test  ban  restrictions  have  prevented  extending 
the  data  base.  Nonetheless,  as  the  present  research  indicates,  there  is  much  more  that 
can  be  learned  from  the  existing  atmospheric  test  data  base. 

To  investigate  the  effect  of  differing  size  distributions  on  both  the  magnitude  and 
longevity  of  optical  effects,  three  distributions  were  selected  for  initial  comparisons— 
one  (the  Defense  Land  Fallout  Interpretive  Code  or  DELFIC  nominal  distribution) 
representative  of  ground  sample  data,  the  other  two  (TTAPS,  R-K)  representative  of 
cloud  sample  data.  A  discussion  of  the  characteristics  of  each  size  distribution  fol¬ 
lows. 


C.l.  DELFIC  Nominal  Distribution. 

DELFIC  is  a  full  physics,  main  frame  fallout  code  developed  by  Hillyer  Norment 
for  the  Defense  Nuclear  Agency  (8).  Although  any  size  distribution  may  be  used  as 
input,  a  nominal  or  default  unimodal  lognormal  distribution  is  provided  which  is 
based  upon  close-in  ground  fallout  samples  over  the  area  out  to  the  10  rad/hr 
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contour  from  events  Small  Boy  and  Ess  (33).  Although  his  n(r)  analysts  has  never 
been  documented,  Norment  says  that  the  close-in  size  distributions  from  these  two 
events  were  very  similar.  This  is  somewhat  surprising  since  Small  Boy  was  a  near- 
surface  burst  and  Ess  was  buried.  A  graph  of  the  DELFIC  size  distribution  and  its 
2nd  moment  (proportional  to  surface  area)  and  3rd  moment  (proportional  to  mass) 
appears  In  figure  11.  The  median  radius,  rm,  and  the  logarithmic  slope  are  .204 /j  and 
4  respectively  (rm  will  always  be  used  to  designate  median  radius  for  the  number  dis¬ 
tribution).  The  distinguishing  characteristic  of  this  distribution  is  its  broad  disper¬ 
sion  (large  logarithmic  slope).  The  median  radius  is  quite  similar  to  those  of  the 
TTAPS  or  E-K  distributions.  The  submicron  particle  population  accounts  for  5%  of 
the  aggregate  surface  area  of  the  DELFIC  distribution. 


C.2.  Nathans/TTAPS  Distribution. 


In  contrast  to  the  simple  unimodal  DELFIC  size  distribution,  the  TTAPS  distri¬ 
bution  is  a  hybrid  function  in  which  the  submicron  population  is  fit  to  a  lognormal 
function,  and  the  supermicron  population  is  fit  to  a  power  law  ''tail”: 


rCl/x;  n(r) 


*■>!/*;  «(r ) 


(13) 


(14) 


The  functions  join  at  rt  which  is  chosen  such  that  the  functions  and  their  first 
derivatives  are  continuous  at  the  splice: 

r0  -  rmesxp  [(P-l)/?2]  (15) 


The  function  was  originally  used  by  Nathans  to  fit  his  surface  burst  cloud  data 
(34).  TTAPS  adopted  Nathans  function  using  parameters  typical  of  Nathans’ 
results:  rm«.25/z,  log  slopes 2,  p«4  and  1.06/x.  The  TTAPS  number,  area,  and 
mass  distributions  are  plotted  in  figure  12.  Besides  being  more  difficult  to  manipu¬ 
late,  the  third  moment  of  the  power  law  tail  (if  p<4)  is  not  analytic  unless  a  max¬ 
imum  particle  radius  is  defined.  If  p™4, 


Figure  11.  DEIFIC  Nominal  Distribution 
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Figure  12.  TTAPS  Distribution 
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where  /m  is  the  mass  fraction  below  r0  (see  Appendix  C  for  a  more  general  expres¬ 
sion  for  rmiix).  The  TTAPS  choice  of  parameters  (r0  — 1.08/i  and  /  «. 084)  yields 
rm„— j  1.3  cm.  The  fraction  of  surface  area  on  particles  below  1  /i  (roughly  r0)  is 
given  by: 
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where  CNF  denotes  the  cumulative  normal  function: 
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and  r2  is  the  median  radius  of  the  second  moment  of  the  submicron  population.  A 
derivation  of  the  Nathans’  function  and  its  moments  is  included  in  Appendix  C. 

The  submicron  population  acoounts  for  60%  of  the  total  surface  area  of  the 
TTAPS  distribution.  This  large  submicron  surface  fraction  yields  optical  effects 
which  are  dramatically  different  in  both  magnitude  and  duration  from  those 
predicted  using  the  DELFIC  nominal  distribution  (figure  3). 


C.3.  Ramaawamy-Kiehl  (R-K)  Distribution. 

The  R-K  distribution,  like  the  DELFIC  distribution,  is  a  unimodal  lognormal 
function.  The  R-K  distribution,  however,  is  more  typical  of  particle  populations 
observed  in  air  burst  clouds  (17).  Ramaswamy  uses  a  surface  mode  radius  of  .26/i 
which  translates  into  a  median  radius  for  the  number  distribution  of  .16/i.  He  uses  a 
logarithmic  slope  of  2.  Nathans  has  overlaid  his  results  for  multiple  air  burst  clouds 
and  determined  that  at  a  median  radius  of  around  .07 /i  and  slope  of  2.1  are  reason¬ 
able  averages  (35). 

The  submicron  population  of  the  R-K  distribution  accounts  for  90%  of  the  sur¬ 
face  area,  larger  even  than  the  TTAPS  distribution.  Thus  the  R-K  distribution 
predicts  the  largest  initial  optical  thickness  and  the  slowest  decay  with  time. 
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D.  Comparative  Results. 

To  illustrate  the  relative  magnitude  and  duration  of  optical  effects  among  the 
DELFIC,  TTAPS,  and  R-K  distributions,  OT  vs  time  was  computed  for  1012  KG  of 
dust  distributed  c  ir  a  hemisphere  (estimate  based  on  5000  MT  exchange).  The  ini* 
tial  cloud  vertical  centroid  was  taken  as  25  km  (peak  density  altitude  of  the  TTAPS 
study).  For  a  description  of  the  code  used  to  compute  stratospheric  debris  removal 
see  section  IV A. 3.  The  results  (figure  3)  clearly  demonstrate  the  large  differences  in 
optical  effects  predicted  by  the  three  distributions.  Note  that  the  TTAPS  and  R-K 
distributions  are  removed  at  nearly  identical  rates  (a  residence  half  life  of  roughly  0 
months).  This  is  explained  by  the  large  fraction  of  surface  area  below  1  micron  which 
is  removed  by  the  turbulent  diffusion  process  rather  than  sedimentation.  The  broad 
DELFIC  distribution,  with  its  large  fraction  of  particles  above  1/u,  is  removed 
predominantly  by  sedimentation  which  accounts  for  its  rapid  decrease  with  time  (half 
life  roughly  1  week).  The  differences  in  OT  behavior  between  the  R-K  and  DELFIC 
distributions  are  predominantly  a  result  of  the  factor  of  two  difference  in  slope. 

E.  Summary. 

It  is  obvious  that  the  optical  occlusion  of  sunlight  energy  responsible  for  climate 
effects  is  extremely  sensitive  to  the  modeled  size  distribution  and  that  several  plausi¬ 
ble  3ize  distributions  would  produce  effects  ranging  from  scant  to  severe.  It  is  also 
clear  that  the  submicron  population  governs  both  the  magnitude  and  duration  of  opt¬ 
ical  attenuation.  The  next  chapter  discusses  activity  injection  into  the  atmosphere 
and  the  data  types  available  for  use  in  the  activity  tracing  model. 
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A.  Data  Records. 

The  radioactive  behavior  of  debris  from  nuclear  events  can  be  used  to  trace  the 
rate  of  its  removal  from  the  atmosphere.  Since  sedimentation  is  a  function  of  particle 
size,  it  is  possible  to  derive  particle  size  statistics  from  the  observed  bulk  vertical 
transfer  of  radioactivity  in  the  atmosphere.  Two  basic  types  of  data  are  available: 

Stratospheric  Removal.  The  stratospheric  burden  of  certain  radionuclides 
has  been  recorded  over  long  periods  of  time  (0,  10,  11,  12,  36,  37,  38,  39,  40,  41,  42, 
43,  44,  45,  46).  Best  records  were  kept  for  90 Sr ,  9sZr,  HTO,  liC  and  185  W.  **Sr 
burdens  have  been  recorded  since  1954  and  reflect  the  contributions  of  92  high  yield 
(megaton  class  and  greater)  events  (38). 

Ground  Deposition.  The  downrange  cloud  arrival  time  was  recorded  for  30 
events  at  the  Nevada  Test  Site  (47).  The  yields  of  these  events  ranged  from  .1  to  50 
KT.  From  this  information  it  was  possible  to  determine  activity  grounded  as  a  func¬ 
tion  of  time.  In  addition,  following  the  1958  test  moratorium,  rough  90 Sr  contours 
were  plotted  over  the  continental  United  States  which  reflect  the  contributions  of 
roughly  75  Nevada  shot3  up  until  that  time  (36). 

B.  Merits  and  Drawbacks  of  Activity  Tracing. 

There  are  major  drawbacks  to  using  activity  tracing.  Unlike  microscopic  tech¬ 
niques  where  actual  debris  is  examined  in  detail,  this  is  an  indirect  method.  In 
order  to  extract  size  information  from  activity  data,  it  was  necessary  to  relate 
radioactive  content  to  particle  size  via  specific  activity  behavior.  Nathans  found  that 
specific  activity  could  in  fact  be  correlated  with  of  particle  size  using  a  single  correla¬ 
tion  function  (17:82).  Unfortunately,  specific  activity  data  is  somewhat  sparse,  highly 
variable  and,  in  some  cases,  difGcult  to  interpret.  This  problem  is  the  weak  link  in 
the  activity  tracing  method  and  may  explain  why  this  approach  has  not  been 
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previously  used  for  particle  size  characterization. 

However,  there  are  some  distinct  advantages  to  the  method  particularly  with 
respect  to  the  nuclear  winter  problem.  Nuclear  winter  is  predicated  on  the  average 
behavior  of  debris  from  hundreds  of  bursts,  and  in  particular,  the  bulk  removal  of 
such  material  from  the  stratosphere.  Indeed,  the  nature  of  the  available  activity 
trace  data  enables  bounding  the  n(r)  behavior  governing  depletion  of  the  stratospheric 
debris  from  a  very  large  number  of  events,  both  U.S.  and  foreign.  Activity  tracing 
data  gives  the  bulk  behavior  of  clouds  and  is  not  subject  to  the  sample  variance 
quandry  associated  with  microscopic  techniques.  Activity  tracing  supplemented  by 
comparisons  with  the  results  of  previous  microscopic  analysis  was  used  successfully  in 
the  present  research  to  determine  n(r)  parameter  bounds  averaged  over  many  events. 

C.  Fallout  Formation  and  Removal. 

An  understanding  of  fallout  formation  and  dispersal  processes  is  important  to 
interpreting  activity  data.  When  a  nuclear  weapon  is  detonated,  fragments  of  fission¬ 
able  material,  unfissioned  active  material,  the  bomb  tamper  and  casing  are  vapor¬ 
ized.  This  collection  of  material  is  referred  to  as  bomb  debris.  Soil  debris  engulfed 
by  the  expanding  fireball  is  also  vaporized  (the  vaporization  temperature  of  soil  is 
'  1700"  K).  Other  close-in  soil  debris  is  melted  or  partially  melted  while  more  distant 
soil  particles  within  the  range  of  the  blast  wave  and  subsequent  gust  is  lofted  in  its 
solid  form.  The  height  (or  depth)  of  the  burst,  the  nature  of  the  terrain,  and  the 
weapon’s  yield  determine  the  relative  amounts  of  soil  debris  which  are  vaporized, 
melted,  or  simply  lofted. 

In  the  case  of  a  free  air  burst,  where  the  rising  fireball  does  not  entrain  surface 
dust>  no  soil  debris  is  present  in  the  highly  radioactive  fireball.  The  free  air  burst 
threshold  altitude  is  '700m  for  1  KT  (13:3-6).  The  fireball  cools  to  debris  vaporiza¬ 
tion  temperatures  within  1-5  seconds.  Within  10  minutes  the  cloud  of  vapor  conden¬ 
sate  stabilizes  at  altitudes  high  enough  (5-35  km  depending  upon  yield)  that  virtually 
no  local  fallout  occurs.  Local  fallout  is  defined  as  that  occurring  within  24  hours  of 
detonation.  The  absence  of  local  fallout  in  the  case  of  free  air  bursts  implies  that  the 
condensed  bomb  debris  particles  are  very  fine.  Indeed,  Nathans  (17,18)  measures  dis¬ 
tributions  tightly  clustered  about  a  median  radius  of  '.1m- 
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For  a  surface  burst  or  low  air  burst  of  a  given  yield,  even  though  the  cloud  rises 
to  similar  altitudes  as  for  a  free  air  burst,  soil  particles  present  provide  the  "carrier" 
for  local  fallout.  The  fact  that  roughly  60%  of  the  activity  falls  within  24  hours 
(48:414)  of  a  surface  detonation  is  due  mainly  to  the  presence  of  larger  particles,  and 
to  a  lesser  degree  to  the  toroidal  hydrodynamic  motion  which  brings  down  a  homo¬ 
geneous  mixture  of  particles  close  to  the  burst  point.  Since  sedimentation  velocity  is 
a  function  of  particle  radius,  it  should  be  possible  to  correlate  activity  down  vs  time, 
A(t),  with  A(r)  hnd  thence  infer  n(r)  from  A(r).  The  larger  particles  are  due  to  the 
presence  of  a  large  mass  of  molten  earth  drawn  into  the  cloud  of  hot  fission  product 
vapor  and  mixed  to  a  greater  or  lesser  degree  by  the  toroidal  circulation  within  the 
cloud.  Judging  from  specific  activity  behavior  for  surface  bursts,  the  mass  of  melted 
or  partially  melted  soil  is  at  least  an  order  of  magnitude  greater  than  the  mass  of 
material  directly  vaporized  (40:388). 

During  the  condensation  phase,  volatile  mass  chains  with  vaporization  tempera¬ 
tures  lower  than  soil  (viz.  As,  Se,  Br,  Kr,  Rb,  Mo,  Tc,  Te,  I,  Xe,  Cs,  ref.  50)  tend  to 
coat  the  surface  of  the  previously  condensed  soil  particles  and  refractory  nuclides 
(those  with  high  condensation  temperatures).  The  refractory  nuclides  (the  products 
of  the  mass  chains  not  listed  above),  tend  to  be  volumetrically  mixed  in  fallout  parti¬ 
cles.  Thus,  the  relative  concentration  of  volatile  nuclides  goes  roughly  as  the  square 
of  the  particle  radius  and  the  concentration  of  refractory  nuclides  goes  roughly  as  the 
cube.  The  different  nuclide  condensation  temperatures,  different  thermal  histories  of 
material,  and  mixing  inhomogeneities  within  the  cloud  lead  to  nonuniform  concentra¬ 
tions  of  the  different  nuclides  on  fallout  particles.  This  phenomenon  is  known  as 
"fractionation"  and  selves  to  complicate  the  activity  tracing  method.  A  further  com¬ 
plication,  noted  by  Nathans  in  his  analysis  of  air  burst  debris,  is  that  below  l/i,  the 
specific  activity  in  some  cases  grows  very  rapidly  with  decreasing  radius  (roughly  as 
r~3).  The  present  effort  has  produced  evidence  that  this  "specific  activity  catas¬ 
trophe"  is  probably  due  to  the  condensation  of  unmixed  fission  products  to  form  high 
activity  particles  (section  V.C.3). 
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D.  Fallout  Data  Categories. 


Fallout  data  may  be  grouped  into  three  general  categories:  global  fallout,  inter¬ 
mediate  fallout,  and  local  fallout.  By  definition,  local  fallout  is  deposited  within  24 
hours  of  detonation  and  within  roughly  (depending  upon  prevailing  winds)  500  km  of 
ground  zero.  Intermediate  fallout  is  not  as  precisely  defined.  For  this  research  it  was 
defined  as  falling  between  one  and  several  days  following  the  detonation  and  as  depo¬ 
sited  within  roughly  5000  km  of  ground  zero  (continental  scales).  Global  fallout  is 
present  only  for  large  yield  (>  IMT)  weapons  capable  of  lofting  contaminated  debris 
into  the  stratosphere.  Fine  material  in  the  stratosphere  lingers  for  months  and  circu¬ 
lates  the  globe  many  times  before  being  removed  by  sedimentation  and  turbulent 
diffusion.  The  longevity  of  debris  in  the  troposphere  is  limited  by  washout /rainout 
and  thermal  circulation  to  periods  of  less  than  3  weeks  (residence  time  is  function  of 
altitude).  Because  the  rainout  and  washout  mechanisms  are  not  operative  in  the  stra¬ 
tosphere,  the  lifetime  of  submicron  particles  reaches  6-12  months  above  '12  km. 

Each  fallout  data  category  is  roughly  associated  with  a  certain  range  of  particle 
sizes  (figure  13).  These  associations  become  apparent  by  studying  the  graph  in  figure 
10  for  the  eve  where  particles  start  from  an  altitude  of  10  km  ('100  KT  cloud  verti¬ 
cal  centroid'.  SArithin  24  hours  particles  20/x  and  greater  will  have  been  removed  by 
sedimentation  and  local  fallout  data  will  yield  information  on  the  fraction  of  particu¬ 
lates  in  this  range.  Within  1  week  (intermediate  fallout  period),  particles  up  to  '5 
microns  will  have  sedimented.  Beyond  a  week,  only  particles  less  than  5  microns  will 
remain  aloft.  Thus  global  fallout  behavior  reveals  information  on  the  submicron  size 
fraction.  These  ranges  are  only  approximate  and  depend  upon  the  initial  cloud 
height  and  wind  velocities. 

D.l.  Global  (Stratospheric)  Fallout  Data. 

Several  programs  have  been  instituted  over  the  years  by  the  Department  of 
Defense  and  Department  of  Energy  to  sample  radionuclides  in  the  stratosphere. 
Early  estimates  of  stratospheric  burden  were  based  upon  limited  balloon  sampling 
and  global  ground  sampling.  With  the  advent  of  thermonuclear  weapons  capable  of 
lofting  ^’g^ificant  amounts  of  debris  into  the  stratosphere,  the  advisability  of  contin¬ 
ued  testing  demanded  an  evaluation  of  stress  of  this  global  fallout  on  the  earth’s 
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environment  and  ecology.  In  1954,  the  Joint  Chiefs  of  Staff  instituted  efforts  at  the 
Defense  Atomic  Support  Agency  (now  the  Defense  Nuclear  Agency)  to  quantify  the 
concentration  of  radioactivity  in  the  stratosphere  (12).  This  led  to  the  organization 
of  the  DoD  High  Altitude  Sampling  Program  (HASP).  An  element  of  this  program, 
Project  Stardust,  provided  the  stratospheric  tracer  data  used  in  the  present  research. 
Beginning  in  1957  Stardust  flew  up  to  two  U2  aircraft  sampling  missions  per  week. 
Sampling  occurred  at  altitudes  up  to  60,000  feet  over  the  North  and  South  American 
continents.  A  map  of  the  HASP  flight  profiles  appears  in  figure  14. 

Following  the  1962  test  ban,  the  U.S.  Atomic  Energy  Commission  (later  the 
Department  of  Energy)  assumed  responsibility  for  sampling  and  instituted  its  own 
High  Altitude  Sampling  Program  including  the  use  of  balloons  (Project  Ashcan)  and 
B57  aircraft  (Project  Airstream)  for  sampling.  The  balloon  probes  sampled  up  to  35 
km,  while  the  aircraft  sampling  was  between  altitudes  of  15-20  km  with  flight  paths 
following  the  west  coast  of  North  and  South  American  from  72°  N  to  50°  S. 

Sampling  was  accompanied  by  extensive  radiochemical  analysis.  Nuclides  which 
have  yielded  the  most  useful  results  for  determining  the  removal  of  debris  from  the 
stratosphere  include  ^Sr,  95 Zr ,  HTO  (tritiated  water  vapor),  14 C  (C02),  and  185  W. 
Of  these  nuclides,  °°Sr  has  been  studied  in  the  most  detail  because  of  its  long  half  life 
(28.1  years)  and  biological  hazard.  Spacial  and  temporal  averaging  of  the  data  was 
use'4  in  plotting  the  stratospheric  burden  of  90 Sr  shown  in  figure  15  (5i).  Of  particu¬ 
lar  interest  are  the  massive  injections  which  occurred  in  the  autumn  of  1961  during  a 
high  yield  Soviet  *est  series.  Because  of  the  extreme  altitude  of  the  cloud  (center 
around  35  km)  and  a  debris  residence  time  which  spanned  several  years,  the  erratic 
effect  of  seasonal  variations  in  the  vertical  diffusion  constant  was  reduced.  In  addi¬ 
tion,  for  high  clouds,  sedimentation  was  not  totally  masked  by  diffusion  as  it  appears 
to  be  for  debris  near  the  tropopause  (37).  Since  90 Sr  has  noble  gas  precursors  which 
condense  late  in  the  fallout  formation  process  it  is  not  the  best  choice  for  tracing  fal¬ 
lout  removal.  Refractory  nuclides,  such  as  96Zr  or  185  W,  are  expected  to  be  more 
evenly  mixed  in  the  debris.  In  addition  refractory  nuclides  tend  to  be  mixed 
volumetrically  in  debris  particles  and  are  less  subject  to  concentration  variation  due 
to  fractionation.  Although  modeling  905r  removal  provided  interesting  and  useful 
results,  it  was  deemed  prudent  to  also  model  the  bulk  vertical  transfer  of  a  refractory 
nuclide. 
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Figure  14.  HASP  Flight  Tracks  (OoO) 
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Figure  15.  90Sr  stratospheric  Inventory,  1951-1974 
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Figure  10.  Moan  Distribution  of  Stratospheric  1<aW,  Sspt.-Oct.  1968 
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For  several  reasons  185  W  was  selected  as  the  refractory  nuclide  for  this  study. 
The  nuclide  was  injected  into  the  stratosphere  only  during  a  three  month  period  in 
1958  (mean  injection  date  of  1  July).  The  tracer  was  an  activation  product  unique  to 
several  devices  in  the  Hardtack  Series.  The  devices  producing  185  W  were  detonated 
on  barges  containing  silica  sand  (10:102).  Project  Stardust  carefully  monitored 
tungsten  data  to  determine  debris  transfer  mechanisms  and  removal  rates  in  the  stra¬ 
tosphere.  Both  the  vertical  and  lateral  diffusion  rates  of  actual  tungsten  clouds  were 
determined  based  on  isopleths  constructed  from  aircraft  sampling  data  (figures  16,17). 
The  stratospheric  tungsten  burden  was  charted  for  a  14  month  period  following  injec¬ 
tion  (figure  18).  The  tungsten  was  injected  at  tropical  latitudes  where  the  effective 
vertical  diffusion  constant  is  low.  Nonetheless,  the  stratospheric  Tungsten  burden 
decays  relatively  quickly,  falling  to  half  its  initial  value  in  5  months.  This  indicates 
that  sedimentation  is  an  important  factor  in  the  activity  transfer  process  since  the 
diffusion  half-life  is  9  months  on  average.  Mason  et  al  (37)  have  shown  that  tritiated 
water  vapor  (HTO)  and  KZr  from  midlatitude  Chinese  shots  are  removed  at  almost 
identical  rates  from  the  lower  stratosphere  (half  residence  time  of  9  months)  implying 
diffusion  is  the  dominant  factor  in  that  case.  A  discussion  of  stratospheric  removal 
processes  is  included  in  section  IVA.3. 

D.2.  Intermediate  Fallout  Data. 

Intermediate  fallout  data  was  the  least  profuse  of  any  category.  Data  taken  and 
analyzed  during  the  DoD  High  Altitude  Sampling  Program  (36,52)  was  used.  HASP 
computed  the  rate  of  deposition  of  stratospheric  debris  using  ®°5r  soil  concentrations 
collected  from  a  world-wide  network  of  sampling  stations  in  1959  (53).  In  order  to 
discriminate  the  stratospheric  fallout  detected  by  U.S.  stations  it  was  necessary  to 
subtract  the  Nevada  Test  Site  contribution  (all  Nevada  shots  were  less  than  100  KT 
and  did  not  contribute  significantly  to  the  stratospheric  burden).  Thus,  it  was  neces¬ 
sary  to  construct  isopleths  of  the  90 Sr  deposition  downwind  of  the  Nevada  Test  Site 
through  1958.  Combining  measured  soil  concentrations  with  known  rainfall  data 
over  the  region,  the  contour  plot  shown  in  figure  19  was  constructed  (52).  A  com¬ 
bined  total  of  995  KT  fission  yield  is  estimated  to  have  injected  100  kilocuries  of  90 Sr 
through  1968  (very  nearly  0.1  kilocurie/KT  fission  on  average,  ref.  36).  An  estimated 
41  kilocuries  of  90 St  fell  within  the  the  10  millicurie/mi?  contour.  Below  10 
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Figure  18.  ^®5W  Stratospheric  Inventory,  1958-1959 
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FIGURE  19.  *°8r  in  U.S.  Soils  from  Intsrmediats  Fallout 

Souro*:  RoftmnM  SZ 
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01  STANCE  FROH  GZ,  STATUTE  MILES 

Source:  Reference  47 

Figure  20.  Operation  PLUMBBOB  -  Coulomb  B. 
Off-site  dose  rate  contours  In  r/hr  at  H+l  hour. 
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millicuries/mi2  it  was  not  possible  to  discriminate  NTS  fallout  from  the  world-wide 
fallout  background.  The  uncertainty  in  the  intermediate  *°Sr  total  is  estimated  to  be 
25%  (52). 

D.3.  Local  Fallout  Data. 

For  many  of  the  U.S.  Nevada  and  Pacific  shots,  an  effort  was  made  to  plot 
activity  contours  down  wind  from  ground  zero.  An  unclassified  compendium  of  fal¬ 
lout  data  for  U.S.  atmospheric  tests  including  ground  deposition  contours  (when 
available)  has  been  published  by  the  Defense  Nuclear  Agency  (47).  For  several  events 
an  approximate  cloud  location  vs  time  is  included  with  the  fallout  contours.  An 
example  is  shown  in  figure  20.  Unfortunately,  most  Nevada  tests  were  air  bursts  and 
since  large  amounts  of  ground  debris  were  not  entrained  in  the  fireball,  only  a  very 
small  part  of  the  activity  was  deposited  locally.  Davis  (54)  has  attempted  to  deter¬ 
mine  particle  size  distributions  from  local  activity  contours  with  some  limited  success 
(discussion  section  VF). 

£.  Summary. 

This  chapter  has  surveyed  the  available  data  and  discussed  the  intrinsic  advan¬ 
tages  and  disadvantages  of  the  activity  tracing  method.  A  strong  point  in  favor  of 
activity  tracing  is  that  data  is  available  on  the  rate  of  activity  depletion  from  the 
stratosphere  over  a  period  of  more  than  a  decade.  The  data  reflects  the  behavior  of 
debris  from  close  to  100  high  yield  events.  Since  stratospheric  material  is  the  source 
of  prolonged  sunlight  attenuation,  this  data  is  particularly  germane  to  the  nuclear 
winter  problem. 


IV.  Modeling  Global  and  Intermediate  Fallout. 


This  chapter  describes  the  modeling  used  to  predict  fallout,  including  both  stra¬ 
tospheric  depletion  and  ground  deposition.  The  models  used  an  assumed  particle  size, 
n(r),  as  input.  The  input  was  varied  in  a  systematic  way  to  make  the  calculated 
activity  burden  and  deposition  match  the  data  summarized  in  the  last  chapter. 

Two  separate  models  were  required.  For  stratospheric  data  comparisons,  a 
model  of  tracer  burden  above  the  tropopause  as  a  function  of  time  was  developed. 
To  utilize  intermediate  fallout  data  a  model  of  ground  deposition  vs  time  was 
required.  The  basic  ingredient  of  each  model  was  an  algorithm  to  compute  cloud  sed¬ 
imentation  and  diffusion  in  one  dimension. 

Admittedly  the  atmosphere  is  not  one  dimensional  and  thus  any  1-D  parameteri¬ 
zation  will  have  limited  accuracy.  However,  in  order  to  facilitate  multiple  variations 
in  n(r)  parameters,  it  was  necessary  to  simplify  the  physics  of  the  atmospheric 
transfer  model.  There  is  precedent  for  1-D  modeling.  One  dimensional  codes  have 
been  used  extensively  in  estimating  the  environmental  impact  of  atmospheric  pollu¬ 
tants.  Ilauer  has  used  1-D  models  to  track  the  transfer  of  nuclear  debris  from  the 
stratosphere  with  reasonable  success  (38,39).  The  TTAJPS  results  were  based  on  a  1- 
D  model  of  the  atmosphere.  Detailed  descriptions  of  the  debris  removal  models  used 

in  this  study  follow.  , 

/ 

A.  Stratospheric  Transfer  Modeling. 

A.l.  Properties  of  the  Stratosphere. 

Megaton  class  nuclear  detonations  inject  large  portions  of  the  resulting  radioac¬ 
tive  clouds  into  the  stratosphere.  Particle  removal  from  the  stratosphere  is  dom¬ 
inated  by  turbulent  diffusion  and  sedimentation  (38,39).  The  sedimentation  process 
was  modeled  using  Stokes  or  Davics-McDonald  fall  mechanics  depending  upon  the 
particle  size  and  altitude  (discussed  in  appendix  A). 
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The  turbulent  diffusion  process  is  not  yet  totally  understood  but  is  strongly 
inllucnccd  by  the  characteristics  of  the  tropopause  (figure  21).  Because  of  its  stable 
temperature  profile,  there  is  very  little  vertical  convective  motion  in  the  stratosphere 
itself.  This  is  especially  true  at  tropical  latitudes  where  clouds  from  nuclear  explo¬ 
sions  can  remain  intact  for  several  passes  around  the  globe  (47:445)  with  less  than  a 
two  mile  change  in  altitude.  Figure  16  provides  evidence  of  this  behavior.  However, 
there  is  considerable  lateral  transfer  in  the  stratosphere  which  is  influenced  by  the 
gaps  in  the  tropopause  in  each  temperate  zone.  The  gap  regions  are  extremely  tur¬ 
bulent  and  it  is  believed  that  a  considerable  transfer  of  air  between  the  stratosphere 
and  troposphere  occurs  there  (confirmed  by  higher  ground  activity  concentrations  at 
latitudes  underneath  the  gap  region,  ref.  36).  The  gap  region  migrates  north  in  the 
summer  and  south  in  the  winter  creating  a  "peeling"  effect  which  is  strongest  in 
spring  and  autumn  months  as  evidenced  by  more  rapid  global  fallout  during  these 
periods  (37,38,39,40,41,42).  in  addition,  the  tropopause  rises  and  falls  with  the 
ground  temperature  which  also  adds  to  the  seasonal  variation  in  transfer  rates. 

Although  the  vertical  diffusion  process  is  stronger  in  the  polar  region  (because  of 
the  gaps)  than  in  the  tropics,  it  is  reasonable  to  define  an  average  vertical  diffusion 
constant  since  nuclear  injections  spread  quite  rapidly  in  the  lateral  direction.  Lateral 
mixing  coefficients  are  of  the  order  of  10s  m2sec-1  (larger  in  the  polar  regions)  which 
yields  a  meridional  spread  of  40°  in  6  months  time  (9).  Thus,  because  of  their  large 
lateral  span,  clouds  experience  an  effective  diffusivity  which  is  in  effect  averaged  over 
latitude.  Bauer  (38)  has  plotted  the  eddy  diffusivity  profiles  used  in  the  one  dimen¬ 
sional  models  of  several  researchers  (figure  22).  For  this  research,  a  seasonal  average 
eddy  diffusivity  (Kz)  of  0.5  m2  sec-1  was  used  with  a  mean  tropopause  height  of  12 
km.  This  value  of  Kz  predicts  a  diffusivity  half  life  for  submicron  particles  in  the 
lower  stratosphere  (12-25  km)  of  about  9  months  which  agrees  with  the  measured 
values  reported  by  Mason  for  and  I1TO  (37). 

Naturally  occurring  stratospheric  aerosols  are  extremely  tenuous  (the  latent  vert¬ 
ical  optical  depth  at  visible  wavelengths  is  only  .005).  It  is  not  expected  that  their 
presence  significantly  affected  the  formation  or  size  distribution  of  nuclear  debris  par- 
ticulales  during  the  condensation  phase.  Thus  it  is  reasonable  to  assume  that 
radioactive  tracers  were  attached  to  weapon  debris  rather  than  ambient  particulates. 
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Figure  21.  Hypothetical  Model  of  Stratospheric  Mixing  and  Transfer 


A.2.  Cloud  Injection. 


A.2.&.  Injection  Altitude. 


Since  the  model  computes  cloud  behavior  over  &  period  of  weeks  to  months  and 
since  clouds  stabilize  in  less  than  10  minutes,  it  was  reasonable  to  omit  cloud  rise 
physics  and  to  initialize  the  computation  at  cloud  stabilization  time.  At  stabilization 
the  cloud  was  given  a  Gaussian  profile  in  altitude  per  Telegadas  and  Bauer  (38:5-2): 


/w)-/. 


(10) 


where 


/(*>*) 


burden  of  tracer  at  z  and  t 
ambient  air  mass  at  z 


M2) 


(20) 


and  ot  =  2.15  km  for  stratospheric  injections.  The  computations  are  fairly  insensi¬ 
tive  to  the  choice  of  initial  dispersion  (38:5-3)  since  the  total  burden  of  tracer  is  being 
computed  over  months.  Results  are  much  more  sensitive  to  the  injection  height  and 
K,. 


There  are  several  empirical  models  for  cloud  injection  height  vs  yield.  The 
Foley-Rudermaa  model  used  in  the  TTAPS  study  fits  a  power  law  to  visible  cloud 
extremities  observed  for  U.S.  tests  (figure  23).  Fits  for  cloud  top  and  bottom  were 
developed  separately  (43): 


CT  -  21.64  Y 2  km 
CB  -  13.41  Y 2  km 


(21) 

(22) 


One  drawback  of  the  Foley-Ruderman  model  is  that  it  is  based  upon  the  dimensions 
of  the  visible  cloud.  The  cloud  radioactivity  profile  may  not  match  the  visible  cloud 
profile,  and  indeed  Zuni  rocket  sampling  indicates  activity  residing  near  the  lower 
boundary  of  the  visible  cloud  (22:24). 

Seitz  et  al  (44)  estimated  cloud  top  and  bottom  based  upon  the  behavior  of 
radioactive  debris  a  few  days  after  the  1961-62  U.S.  and  Soviet  tests.  Because  data 
was  not  available  much  above  an  altitude  or  24  km,  Seitz  arbitrarily  assigned  cloud 
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Figuro  23.  Nucloor  Cloud  RIm  Hoight  at  a  Function  of  Yiold 


top  altitudes  of  24  km  to  high  yield  events.  Thus  Seitz’  curve  fits  are  systematically 
low,  and  should  be  used  with  caution  (45).  Peterson  (46)  has  fit  the  clouds  from  U.S. 
tropical  tests.  However  these  clouds  are  consistently  lower  than  those  from  similar 
bursts  at  northerly  latitudes.  Chang  et  al  (45)  conclude  that  the  Foley-Ruderman 
model  is  useful  as  an  upper  bound,  while  the  Seitz  model  is  useful  for  a  lower  bound. 
Bauer  (38)  believes  the  Seitz  model  to  be  more  realistic,  but  comes  to  this  conclusion 
based  upon  a  fast  diffusion  calculation  of  .1/j  monosized  particles  in  which  sedlmentar 
tioo  effects  are  negligible. 

The  Foley-Ruderman  fit  was  used  in  the  present  effort  to  set  cloud  injection 
height.  The  fit  is  based  on  actual  observations  of  U.S.  events.  There  is  limited  evi¬ 
dence  that  the  fit  predicts  the  cloud  injection  height  of  the  high  yield  Soviet  shots 
(38:3-5).  The  fact  that  activity  has  been  observed  to  concentrate  low  in  the  visible 
cloud  may  have  been  due  to  sedimentation  prior  to  sampling  time.  In  addition,  since 
the  cloud  was  modeled  as  being  distributed  normally  with  atmospheric  pressure 
(equation  10),  the  modeled  activity  centroid  is  below  the  injection  altitude,  za  (see 
section  IVA.3.b).  Bauer’s  underprediction  of  transport  from  Foley-Ruderman  alti¬ 
tudes  is  overcome  if  a  finite  spread  is  introduced  into  his  particle  size  distribution. 
Indeed,  the  spread  which  best  predicts  bulk  vertical  transfer  is  nearly  identical  to 
microscopically  observed  dispersions  for  air  burst  samples  (see  section  VA). 

A2.b.  Event  Data  Input. 

During  the  1050’s  and  1060’s  07  atmospheric  events  lofted  significant  amounts  of 
905r  into  the  stratosphere.  Yield  information  for  most  of  these  events  is  classified. 
To  keep  the  results  unclassified,  events  which  occurred  in  close  succession  were 
grouped  together.  Bauer  (38)  and  Seitz  (44)  have  published  unclassified  grouping 
schemes.  Bauer's  event  groups  were  used  with  slight  modifications  (tables  I, II).  A 
mean  value  of  .1  kilocurie  ^Sr  was  injected  per  kiloton  of  fission  yield  (36:40, 
50:321).  Bauer’s  estimates  of  fission  fractions  were  used  (38:5-3). 

186  W  was  injected  by  several  events  in  the  late  spring  and  early  summer  of  1958. 
For  classification  reasons,  the  model  treated  these  events  as  a  single  injection  occuring 
on  1  July  with  a  stabilized  cloud  center  at  20  km.  The  stabilized  cloud  altitude  was 
based  on  isopleths  generated  from  U2  measurements  (9). 
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TABLE  Is 


Grouping  of  1950s  Injections  -  nSr 

Mean  Injection  Date 

April  54 

June  56 

Sept  56 

June  57 

June  58 

Oct  58 

Country 

US 

US 

USSR 

USSR 

US/UK 

USSR 

Location 

Tropic 

Tropic 

mid-hi  1st 

raid-hi  lat 

Tropic 

Arctic 

Number  of  Events 

5 

5 

4 

5 

30 

10 

Mean  Yield 

O.CMT 

3.5MT 

1.5MT 

2.4MT 

7MT 

3.6MT 

Estimated  90  Sr  In j. 

2.4MC 

880KC 

300KC 

1.4MC 

1.1MC 

1.1MC 

Foley-Ruderinan  Cloud  Center 

20km 

21km 

20km 

21km 

17km 

23km 

TABLE  II: 


Grouping  of  1960s  Injections  -  90 Sr 

Mean  Injection  Date 

Oct  01 

Nov  01 

Jul  62 

Oct  62 

Oct  62 

Country 

USSR 

USSR 

US 

USSR 

/ 

USSR 

Number  of  Events 

13 

2 

7 

5 

11 

Mean  Yield 

2.8MT 

41.5MT 

5.3MT 

25MT 

0  8MT 

Estimated  80 Sr  luj. 

040KC 

OOOKC 

430KC 

4.GMC 

2.0MC 

Foley-Kuderman  Cloud  Center 

22kin 

38km 

25km 

35kin 

26km 
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A.3.  Tracer  Removal. 


A.3.a.  Sedimentation. 


The  particle  sedimentation  model  is  based  on  that  used  in  the  Defense  Land  Fal¬ 
lout  Interactive  Code  (8).  The  Stokes,  Davies-McDonald,  and  Beard  equations  were 
used  depending  upon  the  value  of  the  Davies  number: 


4/?a(Pi“Pa)^3 
3  rfi 


(23) 


where  p4  is  air  density,  pt  is  particle  density,  g  is  the  acceleration  of  gravity,  d  is 
particle  diameter,  and  17  is  dynamic  viscosity  (see  appendix  A). 


The  atmosphere  was  modeled  according  to  equations  in  the  U.S.  Standard  Atmo¬ 
sphere  publication  (55).  Hopkins  (56)  has  shown  that  the  U.S.  atmosphere  is  also  a 
reasonable  approximation  for  the  atmosphere  over  the  Pacific  test  sites  for  purposes 
of  fallout  modeling. 

Particles  were  treated  as  smooth  spheres  of  density  2600  kg/m3  (16)  settling  in 
still  air.  The  Knudsen-Weber  slip  correction  factor  was  used  (57)  which  gives  terminal 
velocities  slightly  smaller  than  the  commonly  used  Cunningham  slip.  The  Knudsen 
formula  fits  terminal  velocity  data  over  a  wider  range  of  d/X^pp  (MFP  denotes  mean 
free  path). 

Appendix  A  presents  the  fall  velocity  equations  used  and  includes  a  printout  of 
the  fall  velocity  subroutine.  Figures  10  and  24  (fall  velocity  vs  size  and  altitude)  were 
generated  using  this  subroutine. 


A.3.b.  Diffusion. 


The  diffusion  model  used  is  similar  to  Bauer’s  "model  I"  (39).  Diffusion  of  a 
chemically  inert  tracer  of  mixing  ratio  f(z,t)  is  described  by  the  following  diffusion 
equation: 


d_ 

dz 


df(z,t) 

dz 


df  ( z,t ) 
dt 


(24) 


where  na(z)  is  the  ambient  atmosphere’s  number  density,  f(z,t)  is  described  by  eqn.  20 
and  is  the  eddy  diffusion  coefficient.  In  this  model,  was  treated  as  constant 


(.5  in2/sec).  The  boundary  at  the  tropopause  was  treated  as  an  infinite  sink  because 
tropospheric  washout/rainout  is  very  rapid  compared  to  stratospheric  removal 
processes. 

If  the  stratosphere  is  treated  as  isothermal,  i.e., 


"•(*)“  ".exp  [~(x -*«)/*„.!«  ] 


(25) 


where  z0  is  the  cloud  center  altitude,  z seaie  is  the  atmospheric  scale  height,  and  n„  i3 
the  tracer  density  at  z0 ,  then  the  diffusion  equation  (24)  simplifies  to: 


K  d2/^*)  _  K*  dj  jz,t) 


dz 2  z scale  dz  '  dt 

Solving  the  equation  for  a  point  source  injection,  we  find: 

r2 


(26) 


/  (*,0 


1 


v2n 


-exp 
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z sc  ale 
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where 


and 


A  **  z0—z 


at  y/2Kdt 

This  f(z,t)  is  equivalent  to  an  n4(z,t)  of  the  following  form: 


(28) 


(29) 


n,(*,0  =  % 


v2i r  o. 


A 


(30) 


which  may  be  thought  of  as  the  product  of  two  probability  distributions,  the  first 
relating  to  the  diffusion  profile  of  the  tracer,  the  second  relating  to  the  atmospheric 
density  gradient.  KN  is  a  normalization  factor  which  is  a  function  of  at  and  z tcale. 
With  some  furthci  manipulation  it  can  be  shown  that  n,(z,t)  is  also  Gaussian: 


I(ff 

V2^rcxp 


z  — 

2 

Z°  7 

,  ‘•scale  , 

l 

(31) 


The  ambient  atmospheric  number  density  gradient  effectively  displaces  the  cloud 
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center  below  z0  by  the  factor  <Ta/zecalt.  As  cr,  increases  with  time,  the  effective  cloud 
center  moves  downward.  Note  that  n4  has  the  same  dispersion  as  f(z,t). 

A.3.C.  Transport. 

Sedimentation  and  diffusion  were  combined  in  a  multigroup  transport  computa¬ 
tion.  The  cloud  was  separated  into  60  velocity  (or  size)  groups: 

(«V (*#»*#  W)  (32) 

where  vg  is  group  sedimentation  velocity  and  zg  is  the  vertical  group  ceiitroid.  The 
cloud  was  treated  in  effect  as  sixty  rigid  clouds  with  separate  streaming  and  eddy 
diffusion  characteristics.  The  initial  vertical  dispersion  of  each  group  was  taken  as 
the  dispersion  of  the  composite  cloud  (2.15  km).  Because  eddy  diffusion  is  a  bulk 
process,  the  same  diffusion  constant  (Krf)  was  applied  to  each  group.  It  was  assumed 
that  once  the  particles  were  formed  during  the  thermal  processes  immediately  after 
detonation,  the  particle  size  distribution  stayed  reasonably  constant  in  time,  i.e., 
agglomeration  was  not  included.  Turco  et  al  estimate  the  effects  of  agglomeration  on 
dust  burden  to  be  less  than  10%  (1:63).  Yoon  et  al  state  that  no  clear  evidence  of 
nuclear  debris  particle  agglomeration  exists  (13:3-43). 

Streaming  was  governed  by  the  sedimentation  velocity  of  particles  at  altitude  zg. 
Since  the  sedimentation  velocity  is  a  function  of  altitude,  treating  the  clouds  as  rigid 
(all  particles  at  centroid  velocity  regardless  of  altitude)  introduces  some  error  into  the 
streaming  calculation.  Since  dv/dz  increases  with  decreasing  radius  (see  figure  24),  so 
docs  the  sedimentation  error.  An  indication  of  the  magnitude  of  the  error  is  given  in 
figure  25.  The  figure  shows  a  comparison  of  the  stratospheric  burden  predicted  by 
rigid  cloud  descent  vs  the  descent  of  a  muitigroup  cloud  (initially  Gaussian)  of  mono- 
size  particles.  The  clouds  were  given  an  initial  centroid  height  of  25  km  and  a  disper¬ 
sion  of  2.15  km.  For  particles  of  radius  >10 /z,  the  error  is  small  enough  to  be  imper¬ 
ceptible  given  the  large  time  increments  (1-2  weeks)  used  in  the  stratospheric  debris 
transfer  calculations.  In  the  submicron  regime,  the  cloud  dispersion  decreases 
significantly  (the  cloud  "bunches  up")  due  to  the  increasing  fall  velocity  with  altitude. 
This  bundling  effect  is  reflected  in  the  more  sudden  drop  in  burden  for  the  multi- 
group  treatment  than  for  the  rigid  cloud  approximation  (figure  25b).  The 
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sedimentation  half-life  is  roughly  the  same  for  the  two  calculations.  The  early  time 
burden  predicted  by  the  rigid  cloud  approximation  is  accurate  to  within  30%  of  the 
multigroup  treatment.  At  later  times,  the  error  is  mitigated  since  for  small  particles 
turbulent  diffusion  is  the  dominant  removal  mechanism,  especially  as  particles 
approach  .In  and  less  (the  results  In  Chapter  V  show  that  in  admissible  size  distribu¬ 
tions  most  particles  are  submicron  with  a  median  radius  of  ".l/i).  Thus,  the  error 
introduced  by  the  rigid  cloud  approximation  is  negligible  above  10 fx  and  probably  less 
than  a  factor  of  2  across  the  remaining  size  spectrum. 

A-3.d.  Burden  Computation. 

The  transport  was  calculated  using  a  Fortran  subroutine  (subroutine  STRAT- 
FAL  of  Appendix  D)  which  injected  clouds  at  Foley-Ruderman  altitudes  and  com¬ 
puted  subsequent  size  group  diffusion  and  sedimentation.  Figure  26  depicts  the  tran¬ 
sport  of  the  l/i  group  as  computed  by  STRATFAL  for  the  Castle  Bravo  event  (area 
under  curves  normalized  to  1).  The  stratospheric  burden  vs  time,  B(t)  is: 

00 

B(t)  -£  /  Aa(rg,z,t)  dz  (34) 

!  *tnt 

where  the  activity,  A,(r)  is  a  function  of  ny(r)  and  the  specific  activity,  S(r).  A  glo¬ 
bal  mean  tropopause  height  of  12  km  was  used.  The  stratospheric  burden  computed 
in  equatiou  34  can  be  directly  compared  with  atmospheric  test  data  as  depicted  in 
figures  15  and  18. 

A.4.  Specific  Activity  Treatment. 

In  order  to  relate  activity  burden  to  number  size  distribution,  a  functional  rela¬ 
tionship  must  be  determined  relating  the  activity  expected  on  a  particle  to  its  radius. 
Freiling  has  proposed  a  "radial  model"  of  activity  content  in  which  refractory  nuclide 
activity  is  approximately  proportional  to  particle  volume: 

Ar{r)  -  Krti(r)r3  (35) 

and  volatile  nuclide  activity  is  approximately  proportional  to  particle  surface: 

Av{r)  -  Kvn{r)r2  (36) 

where  Kr  and  K*  are  constants  of  proportionality. 
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DIAMETER  -  .lOQE+OI  MICRONS 


Flgun*  20.  Group  Fall  Illustration 


Available  data  on  individual  particle  activity  is  usually  presented  in  the  form  of 
equivalent  fissions  per  gram,  a  quantity  called  specific  activity.  Some  authors  use 
"specific  abundance"  which  is  probably  a  better  term,  although  less  prevalent  in  the 
literature.  An  equivalent  fission  is  defined  as  the  number  of  fissions  that  must  have 
occurred  in  a  device  to  produce  the  amount  of  a  particular  radionuclide  observed  on  a 
particle  (49:382).  There  are  5.02  X  1023  nuclei  in  235  grams  of  and  each  fission 
releases  180  MeV.  Thus  1  ton  (2.8  X  1022  MeV)  of  completely  fissioned  material 
when  completely  mixed  with  1  ton  of  soil  will  yield  particles  with  an  average  specific 
activity  of  1.4  X  1014  fissions/g. 

Observed  particle  activities  vary  between  r2  and  r3  over  the  range  of  particle 
sizes  of  interest  for  cloud  modeling,  conforming  to  Freiling’s  hypothesis.  Some 
anomalous  behavior  has  been  noted  and  explained  for  large  particles  (100  -  1000 a*)  by 
R.  C.  Tompkins  (21).  In  addition,  as  previously  discussed  in  section  IVA.4,  Nathans 
has  observed  a  strong,  non-Freiling  radial  dependence  below  .2a*  in  air  burst  debris. 
Since  this  phenomenon  is  probably  due  to  the  presence  of  unmixed  fission  products,  a 
similar  effect  may  not  be  present  for  surface  burets  where  much  more  mixing  material 
(carrier)  is  available. 

It  was  initially  assumed  that  in  the  cloud  most  of  the  activity  resided  on  parti¬ 
cles  between  .01  and  10  microns  and,  further,  that  the  specific  activity  in  this  size 
range  obeyed  Freiling’s  radial  model.  The  first  assumption  was  justified  by  what  is 
known  about  size  distributions  from  microscopic  analysis.  Since  these  analyses  show 
the  bulk  of  cloud  particle  surface  and  volume  lies  below  10a*.  it  was  expected  that  the 
activity  would  probably  behave  likewise.  The  second  assumption  is  confirmed  from 
measured  specific  activity  data  (see  for  example  figures  27-29).  In  most  cases,  specific 
activity  varied  between  that  expected  for  refractory  and  volatile  species: 


re/  ractory: 

volatile : 


Sr(r)  a 
Sv(r ) 


Ar{r) 

- 3~  C 

Av(r) 


constant 


a 


r 


(37) 

(38) 


Note  from  the  example  plots  that  for  the  same  nuclide,  S(r)  behavior  varies  from 
event  to  event. 
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Figure  27 
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Specific  Activities  of  ^®Sr  and 
In  a  Tower  Burst 
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SPECIFIC  ACTIVITY  (FISSIONS 


W 


Figure  28.  Specific  Activities  of  ^Sr  and  147pm 
in  a  Wilson  Sample 
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SPECIFIC  ACTIVITY  (FISSION$/»g) 


In  the  computation  of  stratospheric  burden,  a  normalized  n(r)  was  converted  to 
A(r)  using: 

A(r)  a  n{r)rl  (39) 

From  this  relationship  Af,  the  fraction  of  activity  in  each  group  was  computed  from: 

rt 

A,  -K„  f  A  (r)r‘  dr  (40) 

'1-1 


where 


K.  - 


exp 


IH'm)  ■+ 


in  the  case  of  lognormal  distributions  or 


Kn 


Vfcr/Se  *ltf,rV2rl~  ^  CNF 


In  r0—ln  rt 

+  ( 

0 

i-p+i  y 

-p+i  _  ..i-p+i 


max 


(41) 

>1 

(42) 


in  the  case  of  Nathans’  distribution.  In  the  expression  above,  p  is  the  exponent  of  the 
power  law  tail  and  rf  is  the  median  radius  of  the  /  th  moment  of  the  log  normal  por¬ 
tion  of  the  distribution: 


(43) 


For  90 Sr ,  l  was  set  to  a  nominal  2.5  based  upon  Nathans’  specific  activity 
graphs.  For  the  highly  refractory  186  W,  l  was  set  to  a  nominal  3.0.  A  sensitivity 
study  was  performed  by  varying  l  between  2  and  3  for  both  nuclides.  This  variation 
altered  the  median  radius  bound  by  less  than  a  factor  of  2  for  a  lognormal  n(r). 


A.5.  Matching  the  Model  to  the  Stratospheric  Data. 

wSr  data  comparisons  were  carried  out  by  main  programs  SHOT50  and 
SHOT60  (1950s  and  1960s  comparisons  were  run  separately)  which  varied  the  n(r) 
parameters  (rm ,  0,  p,  fm )  and  then  called  subroutines  to  set  the  activity  distribution, 
inject  the  clouds  from  the  combined  events,  and  then  compute  activity  removal  with 
time.  Appendix  D  lists  the  optimization  program  for  the  1960s  90 Sr  data  comparis¬ 
ons.  Similar  programs  were  developed  for  1950s  905r  and  1958-59  185  W  data  com¬ 
parisons.  In  all  cases,  optimum  n(r)  parameters  were  determined  by  maximizing  a 
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figure  of  merit  expressed  by: 

x-  |s[s+(i;)-B(l,)]2|  ‘  (44) 

where  B+(tj)  represented  the  measured  stratospheric  burden  data  from  the  last 
chapter  and  B(tj)  represented  the  model  results  (eqn.  34).  The  program  flow  is 
shown  in  figure  30. 

Three  sets  of  optimization  runs  were  performed.  1960s  90  Sr  data  was  optimized 
in  a  separate  run  from  the  1950s  00 Sr .  The  third  set  used  the  ■8S  W  data  from  the 
1958  Hardtack  series.  The  90 Sr  data  was  split  for  two  reasons: 

(1)  1950s  data  was  not  as  uniformly  trustworthy  as 
1960s.  Prior  to  Project  Stardust  (under  which  regular  U2 
sampling  was  conducted  beginning  in  1957),  the  90 Sr  burden 
estimates  were  based  on  ground  samples  and  limited  balloon 
sampling.  Thus,  estimated  stratospheric  burdens  prior  to  1957 
were  subject  to  considerable  uncertainty. 

(2)  The  1950s  high  yield  events  were  mostly  contact 
surface  bursts  (detonated  within  1  meter  of  surface).  In  the 
1960’s,  although  the  exact  burst  heights  were  not  available, 
they  tended  to  be  higher  in  general  to  avoid  local  fallout 
(58).  Thus,  running  separate  optimizations  might  identify 
systematic  differences  in  particle  size  distributions  between 
the  two  eras.  Indeed,  results  indicate  a  larger  size 
distribution  for  the  1950s  events  as  compared  to  the  1960s  but 
the  difference  is  small  enough  that  it  could  be  due  to 
uncertainties  in  the  data  and  model. 


Results  from  stratospheric  n(r)  optimization  runs  are  presented  in  section  V.A. 


B.  Tropospheric  (Local/ In  termed  late)  Fallout  Transfer  Modeling. 


B.1.  Properties  of  the  Troposphere. 

Because  of  the  unstable  temperature  profile  in  the  troposphere,  a  different  verti¬ 
cal  transfer  model  was  needed.  Stronger  vertical  air  currents  exist  and  the  hygros¬ 
copic  nature  of  debris  particles  leads  to  removal  by  washout.  Disregarding  horizontal 
transport,  the  following  equation  governs  particulate  density  in  the  troposphere 
(59:507): 

Ki  —  an  —  bn 2  —  ®  (45) 

Qz *  oz 

where  Kj  is  the  eddy  diffusion  coefficient,  n  is  the  vertical  particle  concentration,  a  is 
the  washout  removal  rate  (a  is  the  reciprocal  of  the  residence  time,  T),  b  'is  the  coagu¬ 
lation  coefficient,  and  w  is  the  sedimentation  velocity.  According  to  Junge,  washout 
is  the  predominant  removal  mechanism  for  naturally  occurring  particulates  (average 
radius  .03/x)  below  5  km  so  that  equation  (45)  may  be  simplified  to: 


By  noting  the  altitude  range  over  which  n  decreases  by  a  factor  of  ten,  Junge  esti¬ 
mates  residence  time  for  particulates  below  5  km  to  be  somewhere  between  1  and  10 
days.  Above  5  km,  removal  by  washout  diminishes  and  the  effective  particulate 
residence  time  increases.  The  TTAPS  study  used  a  tropospheric  residence  time  of  10 
days  for  smoke  however  this  represented  a  global  average.  Since  Junge’s  estimates 
were  based  upon  data  taken  over  the  continental  United  States’  midsection,  the 
present  model  varied  T  between  his  bounds.  Even  though  the  average  initial  cloud 
center  height  was  in  the  vicinity  of  10  km  for  Nevada  shots,  sedimentation  lowers  the 
activity  centroid  to  altitudes  for  which  Junge’s  values  were  considered  reasonable. 
Constant  exponential  washout  was  assumed. 

Tropospheric  sedimentation  was  modeled  using  the  same  algorithms  developed 
for  stratospheric  sedimentation  (Section  IVA.3).  Based  upon  Junge’s  analysis, 
washout  was  assumed  to  dominate  eddy  diffusion  effects. 
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B.2.  Tropospheric  Cloud  Injection  anJ  Removal. 

B.2.a.  90 Sr  Injection. 

The  Foley-Ruderman  model  was  not  needed  for  intermediate  fallout  modeling. 
Since  all  events  considered  were  low  yield  U.S.,  the  measured  cloud  heights  were 
available  (47).  The  initial  cloud  centroid  was  taken  as  the  mean  of  the  reported 
cloud  top  and  bottom. 

Initially,  an  attempt  was  made  to  model  injection  and  removal  event  by  event. 
This  approach  proved  extremely  unwieldy  since  it  was  determined  that  75  bursts  con¬ 
tributed  to  offsite  90  Sr  between  1951  and  1958  (Table  III).  Although  wind  data  over 
the  test  site  was  readily  available  at  shot  time,  properly  modeling  individual  burets 
would  have  required  including  wind  and  rain  data  over  the  central  United  States  for 
3  days  following  each  event,  a  grueling  proposition.  Instead,  the  average  behavior  of 
the  75  bursts  was  determined  and  a  composite  injection  of  90 Sr  was  modeled.  Table 
HI  lists  the  contributing  events  and  their  injection  parameters  (based  on  information 
in  reference  47). 

The  75  events  had  a  combined  total  yield  of  995  KT  and  lofted  roughly  99.5 
kilocuries  of  MSr.  Yield  weighted  average  injection  parameters  were  as  follows: 

Y  =  13  KT 

Cloud  center  »•=  10.4  ±1.2  km 

Wind  Speed  -=  70  ±  15  km/hr  (at  cloud  center) 

Wind  Direction  =  246°  ±  24°  (compass) 


B.2.b.  90 Sr  Removal. 


The  model  included  both  sedimentation  (same  model  as  for  stratospheric  tracers) 
and  washout.  The  change  in  the  activity  aloft  as  a  function  of  time  was  described  as 
follows: 


dA 

it 


) total 


—XA  (t)  - 


(  \ 


dA 


(47) 


where  X  is  the  reciprocal  of  the  washout  residency  time,  T.  The  solution  is  of  the 
form: 


TABLE m 

Brente  Contributing  to  09-Site  Fallout  Through  1968 

EVENT 

NAME 

YIELD 

(kt) 

CLOUD 

CENTER  (kfl) 

WIND 

DIRECTION  (•) 

WIND 

SPEED  (mph) 

BURST 

HEIGHT  (ft) 

OPERATION: 

RANGER 

AM* 

1 

13  • 

300 

28 

1060 

Baker-1 

8 

33  * 

290 

38 

1060 

Rmj 

I 

10* 

340 

32 

1060 

Baker-2 

8 

26 

290 

61 

1100 

Fox 

22 

36 

290 

62 

1436 

OPERATION: 

BUSTER- JANGLE 

Baker 

8.5 

27 

60 

24 

1118 

Charlie 

14 

34 

30 

20 

1132 

Dog 

21 

89 

330 

66 

1417 

Eaaj 

31 

43 

340 

62 

1314 

Sugar 

1.2 

13 

210 

61 

3.6 

Unde 

1.2 

9  • 

220 

24 

-17 

OPERATION: 

TUMBLER- SNAPPER 

Able 

1 

13* 

260 

17 

793 

Baker 

1 

13 

340 

9 

1109 

Charlie 

31 

37 

290 

17 

3447 

Dog 

19 

30 

260 

47 

1040 

Eaej 

12 

30  • 

220 

107 

300 

Fox 

11 

37  • 

240 

40 

300 

George 

16 

34  • 

190 

41 

300 

How 

14 

34  • 

160 

29 

300 

OPERATION: 

UPSHOT-KNOTHOLE 

Annie 

IS 

36 

280 

61 

300 

Nancy 

24 

34 

210 

60 

300 

Ruth 

.2 

12 

300 

17 

306 

Dixie 

11 

39 

290 

138 

6022 

Ra, 

.2 

10 

340 

20 

100 

Badger 

23 

30 

310 

63 

300 

Simon 

43 

38 

270 

49 

300 

Encore 

27 

36 

240 

196 

2423 

Harry 

32 

36 

290 

72 

300 

*eatimated 

1  laV*  *  Vi V»WWtW Wff 


rJf.v/T  .T^r*T«nyr4 
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TABLE  m  continued 

Event*  Contributing  to  Off-Sit*  FUiout  Through  1068 


EVENT 

YIELD 

CLOUD 

WIND 

WIND 

BURST 

NAME 

(kt) 

CENTER  (kft) 

DIRECTION  (*) 

SPEED  (mph) 

HEIGHT  (ft) 

II  ~ll 

|  OPERATION: 

UPSHOT-KNOTHOLE 

Orable 

16 

20 

220 

02 

624 

CUmmx 

01 

30 

280 

26 

1334 

OPERATION: 

TEAPOT 

Weep 

1 

18 

320 

78 

782 

Moth 

2 

20 

300 

82 

300 

T«U 

r 

24 

280 

33 

300 

Turk 

43 

40 

200 

18 

600 

Hornet 

4 

33 

280 

43 

300 

Bm 

8 

36 

310 

47 

600  . 

Em 

1 

10  • 

340 

20 

-87 

Apple-1 

14 

7 

280 

47 

600 

Wap' 

3 

27  * 

260 

00 

730 

HA 

3 

33  • 

320 

36 

Hi  Alt 

Po»t 

2 

14 

360 

8 

300 

Met 

22 

30 

240 

88 

400 

Apple-2 

20 

43 

210 

20 

600 

Zucchini 

28 

33 

200 

77 

600 

OPERATION: 

PLUMBBOB 

Bolttmu 

12 

28 

ISO 

23 

600 

Franklin 

.14 

16 

220 

0 

300 

Wilnon 

10 

30 

220 

20 

600 

Priscilin 

37 

34 

280 

16 

700 

Hood 

74 

42 

210 

28 

1600 

Dinblo 

17 

20 

300 

14 

600 

John 

2 

20  • 

220 

22 

Hi  Alt 

Kepler 

10 

24 

180 

18 

600 

Owen* 

0.7 

28 

220 

20 

600 

StokM 

10 

32 

200 

78 

1600 

Shwtn 

17 

24 

300 

7 

600 

Doppler 

11 

31 

230 

43 

1600 

Franklin' 

4.7 

27 

220 

40 

760 

Smoky 

44 

36  • 

300 

36 

700 

OnUleo 

11 

27 

70“ 

12 

500 

■eetimnted 

••extremely  erratic  wind* 
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TABLE  in  continued 

Emu  Contributing  to  09-Sit*  Fallout  Through  1968 


EVENT 

YIELD 

CLOUD 

WIND 

WIND 

BURST 

NAME 

(kt) 

CENTER  (kf».) 

DIRECTION  (•) 

SPEED  (mph) 

HEIGHT  (ft) 

OPERATION: 

PLUMBBOB 

Wb**l«r 

.3 

16 

130 

120 

600 

Coulomb-B 

J 

16  • 

100 

8 

3 

Lop  lac* 

1 

17 

200 

9 

760 

FtMOU 

11 

34 

130 

22 

600 

N*wtoa 

13 

26 

260 

61 

1600 

Wbitoty 

19 

34 

60 

10 

600 

Cbariwton 

13 

26 

100 

48 

1600 

Morg  an 

8 

33 

280 

60 

600 

OPERATION: 

HARDTACK D 

Mora 

3 

14 

360 

23 

1600 

L*a 

1.4 

14 

110 

3 

1600 

Socorro 

8 

23 

220 

19 

1460 

Wran|«ll 

.115 

8 

220 

16 

1600 

Rusbmor* 

.188 

140 

5 

600 

Sanford 

4.9 

19 

230 

32 

1600 

Do  Baca 

2.2 

14 

280 

13 

1600 

Santa  F« 

1.3 

16 

40 

43 

1500 

iT«n|«: 

13.3  kt 

26.2  left 

238' 

39  mpb 

99.5  ft 

*«timat*d 


A{t)  —  A0e 


where  A„  is  the  initial  atmospheric  activity  burden.  A  separate  expression  is  needed 
for  sedimentation  rate  since  this  is  a  function  of  particle  radius  (washout  rate  is  not). 
If  g(t)  is  defined  as  the  fractional  activity  removal  rate  by  sedimentation  such  that: 


where  a  is  normalized  activity  and  A,  as  the  activity  remaining  in  the  atmosphere 
from  which  sedimentation  can  still  progress  (the  "sedimentation  pool"),  then  A(t) 
becomes: 


A(t)-c~*  U0  -/e“  At(<?  )g{t  )  dH 


But  since  the  sedimentation  pool  is  being  exponentially  depleted  by  washout: 


Af  *•  Ag  c 


we  have  finally: 


A(t)  -  e~XtA0  1  -  fg{(  )  it 


This  expression  was  incorporated  into  Fortran  program  INTOPT  (Appendix  E)  as  the 
basis  for  intermediate  fallout  removal  calculations. 


B.3.  Ground  Contour  Computation. 

Equation  52  gives  the  rate  of  tracer  removal.  In  order  to  compare  with  the 
available  data  (which  is  presented  in  the  form  of  ground  contours),  fallout  smearing 
was  modeled.  The  Bridgman-Bigelow  approximation  (29:215)  was  adapted  for  this 
purpose.  Bridgman  and  Bigelow  developed  a  simple  expression  for  ground  activity 
surface  density  (Af)  which  in  terms  of  unit  time  reference  activity  (Aj)  reduces  to: 

,  ,  x  Ai 

Aix.y)- -  (53) 
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where  the  cloud  is  moving  in  the  x  direction,  y  is  the  transverse  direction,  ta  is  the 
cloud  arrival  time,  and  f(y,ta )  is  the  cloud  spacial  distribution  in  the  transverse  direc¬ 
tion  (assumed  to  be  Gaussian).  In  order  to  include  washout,  the  present  model 


replaced  the  Bridgman-Bigelow  g(t)  function  with 


(equation  47).  Other 


otal 


modifications  included  ascribing  a  Gaussian  vertical  profile  to  the  cloud  and  incor¬ 
porating  DELFIC  fall  mechanics  (as  described  in  section  IVA.3.).  A  constant  wind 
velocity  was  assumed  but  the  wind  direction  was  changed  according  to  the  direction 
of  the  activity  hotline  (figure  31). 

Given  Ap(x,y),  contours  were  drawn  using  the  linear  interpolation  routine  con¬ 
tained  in  the  DISSPLAtm  graphics  system  (60).  The  model  was  checked  by  compar¬ 
ing  contours  with  an  example  case  computed  by  Bridgman  and  Bigelow  (29:215-216). 

B.4.  Matching  the  Model  to  the  Intermediate  Fallout  Data. 

For  reasons  discussed  in  section  IV .B.,  rigorous  contour  correlations  between  the 
data  and  the  model  were  not  performed.  Instead,  given  the  tremendous  uncertainties 
inherent  in  the  calculation,  a  simple  comparison  of  activity  grounded  vs  distance  was 
used  in  an  attempt  to  optimize  n(r).  The  activity  contours  were  subdivided  by 
breaking  the  hotline  into  six  segments  as  shown  in  figure  31.  The  activity  deposited 
along  each  hotline  segment  was  estimated  assuming  a  Gaussian  transverse  distribu¬ 
tion  such  that: 


AA,  *»  AP(Ax,V27r  <ry> 


(54) 


where  APi  is  the  peak  surface  activity  density  along  hotline  segment  i,  Ax,  is  the  seg¬ 
ment  length,  and  was  estimated  from: 

a  W/2 


°jr.  %  ^io 


2  In- 


10 


(55) 


Yl0  was  determined  by  bisecting  each  hotline  segment  and  averaging  the  distance 
from  the  hotline  to  the  10  mC/mi2  contour.  While  upper  and  lower  bounds  for  the 
activity  density  were  available  for  segments  2-6,  no  upper  bound  was  available  for 
segment  1.  Possibly  high  values  of  activity  density  due  to  local  fallout  within  the 
60  mC/mi2  contour  would  have  made  a  large  difference  in  the  average  activity 
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TABLE  IV 


Incremental  Intermediate  Fallout  Deposition 

segment 

Ax 

AA 

number 

millicur/mi2 

miles 

miles 

kilocuries 

1 

65 

459 

190 

14.7 

2 

55 

757 

310 

32.3 

3 

45 

252 

284 

8.08 

4 

35 

229 

253 

5.09 

5 

25 

184 

237 

warn 

0 

15 

23 

394 

0.34 
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density  along  segment  1.  Checks  of  local  maps  of  90 Sr  deposition  in  Nevada  and 
Utah  revealed  that  85  mC /mi2  is  a  reasonable  average  (61,62).  Table  IV  gives  the 
deposition  computed  tor  each  segment. 

This  method  allows  for  the  presence  of  activity  outside  the  10  mC/mi2  contour 
(in  contrast  to  roughly  40%  computed  within  the  10  mC  /mi2).  Isotopes,  Inc.  (52) 
constructed  the  intermediate  fallout  map  used  for  the  analysis  and  guessed  that  the 
bulk  of  the  remaining  60%  was  depo  d  locally.  The  present  analysis  supports  the 
opposite  conclusion,  namely  that  the  remaining  activity  was  still  aloft. 

The  intermediate  fallout  n(r)  optimization  proceeded  similarly  to  the  stratos¬ 
pheric  n(r)  optimization  (Section  IV.A.5).  The  only  difference  was  that  the  intermedi¬ 
ate  fallout  n(r)  optimization  was  based  on  ground  activity  fraction  rather  than  the 
stratospheric  fraction.  Optimization  was  performed  by  a  main  program  INTOPT 
which  varied  n(r)  parameters,  called  subroutine  INTFAL  to  compute  tracer  removal, 
and  finally  computed  a  figure  of  merit  as  in  equation  44.  The  results  of  the  optimiza¬ 
tion  are  discussed  in  section  V.B. 


V.  Results. 


The  previous  chapter  described  the  modeling  of  vertical  activity  transfer  in  the 
stratosphere  and  troposphere  and  the  approach  used  to  optimize  n(r).  Chapter  V 
presents  the  results  of  the  optimization  calculations.  Chapter  VI  summarizes  the 
findings  and  presents  recommendations  for  their  implementation. 

A.  Global  (Stratospheric)  Tracer  Removed. 


A.1.  *°Sr  Removal. 


The  unaltered  DELFIC  size  distribution  was  used  as  the  initial  n(r)  input  to  the 
stratospheric  fallout  subroutine,  STRATFAL  (see  listing,  Appendix  D).  The  DELFIC 
distribution  was  obviously  too  heavy  to  explain  the  data  as  is  evident  from  figure 
32.  A  combination  of  the  1950s  and  1960s  results  are  shown.  If  the  activity  were  dis¬ 
tributed  according  the  the  DELFIC  ground  fallout  nominal,  the  stratospheric  burden 
would  have  decayed  in  days  rather  than  months.  The  DELFIC  distribution  is  obvi¬ 
ously  biased  much  too  strongly  toward  the  heavy  end  of  the  particle  size  spectrum. 


It  was  expected  that  the  TTAPS  distribution,  with  its  larger  activity  fraction 
below  ip,  would  better  match  the  data.  This  was  indeed  the  case,  as  shown  in  figure 
33.  Although  some  improvement  was  evident  in  comparison  to  the  DELFIC  nominal, 
even  this  distribution  was  too  heavy  to  explain  the  data  (the  modeled  stratospheric 
burden  again  drops  too  quickly).  The  critical  parameter  in  Nathans’  distribution  is 
/m,  the  mass  fraction  below  *•„.  This  turns  out  to  be  a  rather  complicated  function 
of  Nathacs’  n(r)  parameters: 
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Figure  32.  *Sr  Burden  from  U.S.  and  Foreign  Taata,  OELFIC  Nominal 
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Figure  13.  *Sr  Burdan  from  U.9.  and  Foreign  Tasta,  Nithana'  Distribution 


R0  is  a  function  of  the  power  exponent,  p,  and  varies  between  0.7  and  1.7 fi  for 
3<p<5.  Thus  over  the  range  of  p  bracketing  Nathans’  data,  fm  is  a  good  indica¬ 
tor  of  the  fraction  of  mass  in  the  submicron  range. 

Best  fits  to  the  stratospheric  data  were  obtained  when  all  of  the  mass  was  placed 
in  the  submicron  regime  (/m  — ►  1.0).  For  fixed  rmftT,  fm  is  optimized  for  high  values 
of  p  (p  — *  oo).  For  fixed  p,  f  m  is  maximized  as  — *ra .  In  either  case,  best  fits 

to  the  stratospheric  data  were  obtained  with  fm  »  1  such  that  optimum  rmax  and  p 
values  are  not  physically  reasonable.  As  an  example,  figure  34  depicts  the  datarmodel 
comparison  with  100%  of  the  mass  below  r0 .  The  fit  looks  good,  but  the  size  distri¬ 
bution  is  truncated  at  r0.  Thus,  the  Nathans’  hybrid  function  does  not  adequately 
predict  the  removal  of  90 Sr  from  the  stratosphere  . 

Next,  parameters  for  a  pure  lognormal  distribution  were  optimized  by  varying 
rm  from  .01  to  10 fx  and  the  slope  from  1.01  to  4.  A  Gauss-Newton  optimization  was 
used  at  first  but  later  abandoned  because  the  data  did  not  yield  a  unique  solution  (see 
figures  36,  36).  It  was  more  expedient  to  vary  the  parameters  continuously  and  gen¬ 
erate  a  three  dimensional  plot  of  the  figure  of  merit,  x>  vs  Tm  and  slope.  The  figure 
of  merit  was  defined  as: 

Xt'nA  -  js  [s-(l,)  -  '  (57) 

where  B+(<;)  represents  the  measured  stratospheric  burden,  and  B(t)  the  model 
results. 

Separate  parameter  variations  were  performed  for  1950s  and  1960s  data.  A  total 
of  460  cases  were  run  for  each  era.  Figures  35  and  36  plot  x  vs  rm  and  slope  for 
1960s  wSr  removal  and  1960s  wSr  removal  respectively.  Note  the  lack  of  a  unique 
solution  (see  figures  36,  36).  Rather,  for  any  slope,  an  optimum  rm  is  indicated.  If 
the  slope  is  ~2  (Nathans’  value)  the  median  radius  is  in  the  vicinity  of  .In  for  each 
optimization,  somewhat  lower  than  the  rm  indicated  by  Nathans’  surface  burst  sam¬ 
ple  analysis.  The  figure  of  merit  improves  slightly  as  rm  decreases.  Thus  the  distri¬ 
bution  which  best  predicts  the  removal  of  ^Sr  from  the  stratosphere  is  biased  more 
toward  smaller  particles  than  evea  the  TTAPS  distribution  and  can  be  modeled  as  a 
simple  lognormal  function  rather  than  a  hybrid  splice.  Figure  37  plots  model  results 
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Figure  34.  *Sr  Burden  from  U.S.  and  Foreign  Teats,  fm  *  1  (rmax 
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Figure  36.  N(r)  Parameter  Variation.  1960a  *°Sr 
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vs  stratospheric  data  using  rm  —  .Ifx  and  /?  =  ln(2).  The  fit  is  quite  reasonable  given 
the  large  uncertainties  in  the  stratospheric  measurements.  The  disparity  between 
model  and  data  around  1960-61  is  possibly  explained  by  the  arrival  of  debris  from 
two  high  yield  exoatmospheric  bursts  (250  km),  Teak  and  Orange,  which  occurred  in 
late  1958  (52:82).  These  bursts  were  not  included  in  the  burden  calculations. 

It  might  be  argued  that  the  bias  toward  a  smaller  particle  distribution  was  a 
result  of  the  volatility  of  90 Sr  precursors  (63:3).  The  mass  chain  condenses  late 
(roughly  5  minutes  after  detonation  for  megaton  class  weapons)  such  that  the  gaseous 
precursors,  not  heavily  influenced  by  gravity,  may  have  risen  higher  than  early  con¬ 
densing  chains.  The  rare  gas  precursors  would  have  inhibited  90  Sr  scavenging  by 
local  fallout.  At  the  time  of  condensation,  sedimentation  would  have  lowered  the 
concentration  of  the  heaviest  particles  to  some  degree,  such  that  condensation  nuclei 
may  have  been  generally  smaller  for  90 Sr .  The  best  test  of  this  possible  explanation 
for  the  smallness  of  the  particulates  carrying  the  90 Sr  was  to  compare  the  behavior  of 
a  refractory  nuclide.  Fortunately,  data  on  the  stratospheric,  removal  of  185  W,  with 
one  of  the  highest  condensation  temperatures  on  the  periodic  chart  (therefore  highly 
refractory),  were  available. 

A.2.  185  W  Removal. 

The  sequence  of  optimization  trials  for  185  IF  was  identical  to  that  used  for  90 Sr . 
As  in  the  case  of  90 Sr,  the  nominal  DELFIC  distribution  predicted  a  stratospheric 
removal  rate  which  was  much  too  rapid  (figure  38).  But  again,  surprisingly,  the 
TTAPS  distribution  predicted  a  faster  removal  than  indicated  by  stratospheric  sam¬ 
pling  data  (figure  39).  Indeed,  the  behavior  of  185  W  appeared  to  be  very  similar  to 
that  for  "Sr,  such  that  the  TTAPS  function  again  did  not  give  satisfactory  results 
except  for  non-physical  values  of  n(r)  parameters  (see  discussion  in  section  VA.l). 

As  with  "Sr,  parameters  were  optimized  for  a  lognormal  function  based  on  a 
figure  of  merit  as  expressed  in  equation  57.  The  results  were  not  extremely  different 
from  "Sr .  .Apparently,  the  two  nuclides  resided  on  similar  debris  particle  distribu¬ 
tions.  A  reasonable  datarmodel  comparison  for  188  W  is  shown  in  figure  40  for 
rm  —  ,2/i  and /? -• /n(2).  Figure  41  plots  best  fit  rm  and  logarithmic  slope  values 
from  "Sr  and  185  IF  optimization  calculations.  Notice  that  the  185  W  results  are 
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Figure  40.  185^  from  Hardtack  Series, 

Reasonable  Fit 
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Figure  41.  Stratospheric  Strontium  and  Tungsten  Removal 
Best  Fit  Median  Radius  and  Slope 


straddled  by  the  “Sr  results,  evidence  that  the  tracers  reside  on  similar  particle 
populations  despite  diriering  mass  chains  and  physical  properties.  These  curves  by 
themselves  did  not  isolate  a  unique  size  distribution.  In  order  to  limit  the  range  of 
admissible  parameter  values,  another  experimental  observation  was  needed. 


In  his  experience  with  test  sample  data,  Nathans  noted  that  he  rarely  if  ever  saw 
nuclear  debris  particles  of  radius  less  than  .01 /i  (64).  Since  he  based  his  size  distri¬ 
bution  statistics  on  samples  of  typically  1000  particles,  the  probability  of  the  presence 
of  particles  <  .01  is  at  most  .001.  By  plotting  the  locus  of  rm  vs  slope  which  satisfy 
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it  was  possible  to  bound  the  admissible  values  of  rm  and  slope  (figure  42  shaded 
area).  Since  for  each  nuclide  the  figure  of  merit,  improved  with  decreasing  rm ,  actual 
parameters  are  expected  to  lie  closer  to  the  lower  boundary  of  the  admissible  area 
(lower  rm ,  higher  logarithmic  slope). 

The  particles  which  linger  in  the  stratosphere  obey  a  size  distribution  which  may 
be  described  by  a  lognoTnal  function  and  is  smaller  in  general  than  either  the 
DELFIC  or  TTAPS  distributions.  Thus,  stratospheric  particulates  were  initially 
expected  to  be  removed  at  an  even  slower  rate  than  the  nuclear  winter  studies 
predicted  (1,2, 3, 4, 5).  However,  this  conclusion  turned  out  to  be  premature,  and  was 
contradicted  by  the  intermediate  and  local  fallout  results. 


An  important  clue  to  determining  the  true  size  distribution  parameters  was  the 
fact  that,  for  both  90 Sr  and  185  W,  the  size  distribution  which  predicted  their  long 
term  behavior  was  quite  similar  to  size  distributions  which  Nathans  determined  for 
air  bursts.  This  was  despite  the  fact  that  all  bursts  producing  185  W  were  contact 
surface  bursts,  and  that  the  1950s  events  producing  significant  stratospheric  90Sr 
were  predominantly  surface  or  near-surface  bursts. 


B.  Intermediate  Fallout. 

As  described  in  section  IV.B.  intermediate  fallout  data  consisted  of  90 Sr  activity 
contours  over  the  continental  U.S.  The  initial  comparison  calculation  used  the  n(r) 
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bounds  determined  from  the  analysis  of  stratospheric  90  Sr  removal  (nominally  taken 
as  rm  =  .1  fi,  B  =  ln(2)).  The  results  were  quite  puzzling.  Virtually  no  intermediate 
fallout  due  to  sedimentation  was  predicted  if  the  *°Sr  was  carried  by  particles  obey¬ 
ing  a  stratospheric  size  distribution.  Washout  appeared  to  be  the  only  removal 
mechanism  (figure  42).  A  reasonable  comparison  between  predicted  and  actual  con¬ 
tours  was  obtained  by  turning  off  the  sedimentation  effect  entirely  and  using  a 
washout  e-folding  time  of  48  hours  (figure  43). 

Because  of  the  tremendous  uncertainties  inherent  in  the  computation  of  inter¬ 
mediate  fallout  (combining  75  bursts,  assuming  constant  wind  along  the  hot  line, 
assuming  exponential  washout  which  may  completely  mask  sedimentation  effects),  it 
was  not  possible  to  predict  size  distribution  parameters  from  intermediate  fallout 
data  with  any  confidence.  However  one  set  of  optimization  calculations  was 
made  using  plausible  atmospheric  parameters  (constant  wind  speed  of  40  mph, 
washout  e-folding  of  72  hours)  to  determine  if  the  size  distribution  could  be  described 
with  a  lognormal  function.  Indeed  it  could,  as  shown  in  figure  44.  A  median,  radius 
of  .27 n  and  slope  of  3.1  optimized  intermediate  90 Sr  removal  for  the  assumed  atmos¬ 
pheric  behavior. 

As  a  point  of  interest,  the  model  was  exercised  using  the  DELFIC  nominal  distri¬ 
bution  and  the  TTAPS  distribution,  with  results  shown  in  figures  45  and  46.  The 
DELFIC  distribution  obviously  overpredicted  grounded  activity.  Indeed,  the  DELFIC 
distribution  by  itself  came  close  to  predicting  observed  contours.  However,  it  is 
extremely  unlikely  that  the  true  distribution  resembled  the  DELFIC  n(r)  since  most 
of  the  75  bursts  contributing  to  off-site  activity  were  air  bursts  and  the  nominal 
DELFIC  distribution  is  based  on  contact  surface  and  buried  burst  particles.  The 
TTAPS  distribution  gave  a  reasonable  fit  but  this  was  considered  fortuitous  given  the 
uncertainties  in  the  model.  A  rigorous  optimization  of  the  TTAPS  distribution  was 
considered  to  be  unfruitful. 

The  most  useful  information  from  the  intermediate  fallout  analysis  was  that  the 
stratospheric  distribution,  based  largely  upon  data  from  surface  and  near-surface 
burst  injections,  predicted  virtually  no  fallout  by  sedimentation.  This  result  v/as  puz¬ 
zling  since  significant  local  fallout  was  observed  on  clear  days  (no  washout)  for  all 
surface  and  near-surface  U.S.  events  (47).  The  intermediate  fallout  results  implied 
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Figure  43.  Composite  Event  Fallout  Contours 
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Figure  44.  Activity  Grounded  vs.  Distance, 
Best  Fit  n(r) 
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Figure  45.  Activity  Grounded  vs.  Distance 
DELFIC  Nominal. 
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Figure  46.  Activity  Grounded  vs.  Distance, 
Nathans'  Distribution 
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that  the  stratospheric  n(r)  did  not  represent  the  complete  description  of  the  actual 
particle  size  distribution  for  surface  bursts. 

C.  Resolving  the  Global  vs  Local  Fallout  n(r)  Dichotomy. 

The  stratospheric  and  intermediate  fallout  results  were  paradoxical.  The  stratos¬ 
pheric  tracers  resided  on  particles  lofted  predominantly  by  surface  and  near-surface 
bursts.  Yet  the  stratospheric  activity  resided  on  particle  distributions  more  represen¬ 
tative  of  air  bursts.  The  intermediate  fallout  model  predicted  no  activity  deposited 
via  sedimentation  when  upper  bound  stratospheric  size  distributions  were  used  thus 
contradicting  observed  local  fallout  particle  size  distributions  from  surface  bursts 
(DELFIC  being  one  example).  These  local  size  distributions  arc  removed  primarily  by 
sedimentation  as  predicted  by  Bridgman-Bigelow  (29)  and  others.  Additional  infor¬ 
mation  was  sought  to  explain  coexistence  of  these  different  size  distributions  each  of 
which  explained  one  portion  (global  vs  local)  of  the  observed  fallout  but  not  the 
other. 


C.X.  Cloud  Samples  vs  Ground  Samples. 

The  differences  observed  between  particles  from  cloud  and  ground  samples  pro¬ 
vided  an  important  clue  to  resolving  the  size  distribution  dichotomy.  Tompkins  et  al 
have  published  a  comparative  study  (49:381-100)  of  cloud  and  ground  samples  from 
four  events  in  which  they  superimposed  plots  of  specific  activity  (S)  vs  particle  size 
Tor  cloud  and  closc-in  fallout  samples.  Figures  47-50  show  these  superpositions  for 
1,1  Cc,  "Mo,  M0Ba,  and  /,5Ca  respectively.  A  careful  examination,  of  these  plots 

/  f  t 

reveals  that,  invariably,  S(r)  is  lower  for  ground  samples  than  for  cloud/distant  sam¬ 
ples  by  a  significant  margin  (roughly  a  factor  of  5  for  refractory  chains  when  aver¬ 
aged  over  radius;  less  pronounced  for  volatile  chains).  Thus  it  appears  that  for  the 
same  event  specific  activity  is  bivariate,  and  fundamental  differences  exist  in  the  phy¬ 
sical  properties  of  particles  deposited  locally  vs  particles  remaining  in  the  cloud  or 
deposited  al  some  distant  location. 

'Phe  existence  of  these  fundamental  differences  is  further  supported  by  Tompkins’ 
separation  of  ground  samples  from  the  5  megaton  barge-mounted  Tevva  event  (sam¬ 
ples  were  from  a  location  17.5  km  upwind  and  probably  represent  debris  which  feil 


91 


SPECIFIC 


GEOMETRIC  MEAN  PARTICLE  DIAMETER  H 

•  TOwa  Spharaa 

A  Zunl  Distant  Fallout 
▼  Zunl  Cloaa-ln  Fallout 
0  Bravo  Fallout 
O  Zunl  Cloud 
□  Tawa  Irtagulaim 
7  Lacroaaa  Fallout 

RKRRINTCD  WITH  PERMISSION  FROM  ADVANCES  IN  CHtMtSmY  SCRIIS: 
HAOIOMUCVOeS  IN  THE  VtVWONMtNT.  COPYRIGHT  1*70  AMERICAN  CHEMICAL  SOCIETY 
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Figure  48,  Specific  Activity  of  **Mo  in  Fallout  and  147Pm 
in  Cloud  varaua  Particle  Size 
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Figure  49.  Specific  Activity  of  ^Ba  in  Fallout  and  “Sr 
in  Cloud  varaua  Particia  Siza 
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within  1  hour  of  detonation)  into  spherical  and  irregular  fractions.  The  specific 
activity  of  these  fractions  are  overlaid  in  figure  48.  It  is  interesting  that  the  higher 
specific  activity  of  the  spherical  particles  (over  their  admittedly  limited  size  range) 
more  closely  matches  that  of  the  cloud  particles  while  the  lower  specific  activity  of  the 
irregular  particles  more  closely  matches  that  of  the  ground  samples.  Thus  even 
though  Tompkins  did  not  separate  cloud  samples  into  spherical  and  irregular  frac¬ 
tions,  one  might  expect  that  cloud  samples  would  be  richer  in  spherical  particles  and 
close-in  samples  would  be  richer  in  irregular  particles.  Tompkins  concluded  that  the 
high  specific  activity  of  the  spherical  particles  was  the  result  of  a  more  intense  ther¬ 
mal  history  such  that  refractory  chains  were  incorporated  in  the  spherical  particles 
during  a  hot  early  phase  which  the  irregular  particles  never  experienced  (49:395). 
Apparently  the  spheres  originated  from  the  condensation  of  material  vaporized  in  the 
fireball  at  a  time  when  the  fission  product  concentration  was  very  high.  The  irregu¬ 
lar  particles  did  not  reach  fireball  temperatures  and  probably  were  unmelted  or  par¬ 
tially  melted  soil  particles  whose  material  was  not  homogeneously  mixed  with  fission 
products,  hence  they  exhibited  a  lower  specific  activity. 

C.2.  Spherical  and  Irregular  Particles  in  Cloud  Samples. 

Nathans  confirmed  the  existence  of  spherical  and  irregular  particles  in  debris 
clouds  by  separating  cloud  samples  into  spherical  and  irregular  particles  on  two  sili¬ 
cate  surface  bursts  of  yield  .5  KT  and  10  KT  (16:307 ;64).  He  found  that  his  smaller 
size  fractions  were  richer  in  spherical  particles,  and  that  the  bulk  of  radioactivity  was 
carried  by  the  spheres  (as  a  consequence  of  their  higher  specific  activity).  For  the  10 
KT  event,  he  performed  independent  size  statistics  for  spherical  and  irregular  parti¬ 
cles  and  found  that  they  obeyed  different  lognormal  distributions  (figure  51).  The 
spherical  particles  tended  to  be  smaller  with  median  radius  of  0.1  n  and  a  logarithmic 
slope  of  2.  The  irregular  particles  tended  to  be  larger  and  more  disperse  with  median 
radius  of  .17 /i  and  a  logarithmic  slope  of  3.  For  both  events  he  found  that  spheres 
outnumbered  irregulars  in  cloud  samples  approximately  2.3:1.  The  larger  irregular 
particle  ensemble  would  sediment  more  rapidly  than  the  smaller  more  narrow  spheri¬ 
cal  particle  distribution.  It  is  interesting  that  the  size  distribution  of  spherical  parti¬ 
cles  is  typical  of  the  size  distributions  Nathans  measured  for  air  bursts  (17).  This 
fact,  taken  together  with  the  optimum  distribution  derived  for  stratospheric  debris, 
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implies  that  surface  burst  particle  distributions  contain  a  component  which  is  similar 
if  not  identical  to  the  size  distribution  for  free  air  bursts  (in  which  little  or  no  soil  is 
mixed  with  the  weapon  debris  and  all  debris  reaches  fireball  temperature). 

C.3.  Specific  Activity  Catastrophe. 

From  free  air  burst  data  it  is  apparent  that  debris  which  reaches  fireball  tem¬ 
perature  tends,  through  hydrodynamic  processes,  toward  a  characteristic  size  distribu¬ 
tion  with  rm  «s  l/i,  0  «  ln( 2).  This  behavior  is  corroborated  by  the  size  distribution 
of  what  appears  to  be  unmixed  fission  products  in  air  burst  debris.  Nathans  observed 
that  below  l/i  the  specific  activity  of  particles  increased  rapidly  with  decreasing 
radius  (64)  in  a  manner  uncharacteristic  of  the  Freiling  radial  behavior  (in  which 
activity  increases  as  the  2nd  or  3rd  moment  of  the  number  size  distribution,  refs.  66, 
67).  An  example  of  this  behavior  for  Promethium  is  shown  in  figure  52.  Nathans 
described  the  effect,  which  leads  to  extremely  high  specific  activities  at  low  radius,  as 
a  "catastrophe"  which  he  was  unable  to  explain.  For  the  several  airburets  which  he 
studied,  the  specific  activity  for  submicron  particles  was  observed  to  vary  between 
r“25  and  r“35  for  refractory  chains. 

The  effect  was  of  concern  for  the  present  research  because  the  fraction  of  total 
activity  carried  by  submicron  particles  could  be  dramatically  affected  by  the  catas¬ 
trophe.  The  effect  can  be  explained  if  at  the  time  of  detonation  a  totally  random 
fission  product  dispersion  occurs  such  that  nuclide  concentrations  on  individual  parti¬ 
cles  are  roughly  independent  of  radius.  Particle  specific  activity  would  then  be  scat¬ 
tered  about  r-3,  as  observed.  In  one  study  (64)  Nathans  overlaid  95/r  specific 
activity  data  for  a  range  of  yields.  A  linear  fit  to  his  points  in  log  space  reveals  that 
the  submicron  specific  activity  obeys  the  following  relation  roughly  independent  of 
yield: 

^  16  X  1017  /  i 8810116  ,  ^ 

S(r )  - — - - -  (r  in  mtcrone )  (59) 

j-2  5  gram 

It  i3  hypothesized  that  the  submicron  particles  are  pure  fission  products,  and  that 
they  therefore  obey  the  characteristic  size  distribution  for  high  temperature  debris.  If 
the  fission  product  dispersion  is  totally  random,  it  is  useful  to  compute  an  average 
activity  per  gram.  Tbb  average,  if  indeed  the  particles  are  composed  of  pure 
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Figure  52.  Specific  Activity  Catastrophe  Example 


fissioned  material,  should  equal  the  total  activity  available  (1.4  X 1026 
fissions/megaton)  divided  by  the  mass  of  totally  fissioned  bomb  material  (5.7  X  104 
g/megaton).  Using  the  number  size  distribution  characteristic  of  stratospheric  surface 
burst  debris  and  air  burst  debris,  the  average  specific  activity  is  in  fact  within  10%  of 
this  quotient,  i.e.: 


<S(r)> 


OO 


/  s(r)n(r)dr 


1.4  X  1026  fissions 
5.7  X  104  (Tam 


(60) 


where, 

*(r)~  2)rCXP 

The  equality  expressed  by  equation  80  lends  further  credence  to  the  parameter  values 
believed  to  govern  the  high  temperature  debris  size  distribution. 
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D.  Proposed  Cloud  Distribution. 

The  evidence  discussed  in  section  VA,  V.B,  and  V.C  supports  the  existence  of  a 
bimodal  cloud  particle  distribution  for  surface  bursts: 

n(r)  •  nj(r)  +  «2(r)  (62) 

where  n  l  represents  high  activity  spherical  particles  formed  from  the  condensation  of 
vaporized  bomb  and  surface  materials.  The  second  mode,  n2,  represents  irregularly 
shaped  particles  of  unmelted  or  partially  melted  (low  temperature  history)  surface 
material  with  a  lower  activity  concentration  in  general. 

Such  a  distribution  explains  the  apparent  differences  in  size  distributions  govern¬ 
ing  the  behavior  of  fallout  at  early  and  late  times.  In  the  stratosphere  at  late  times, 
an  nt  type  distribution  governs  tracer  removal  as  demonstrated  by  the  90  St  and 
186  W  optimization  results  of  section  V.A.  The  early  cloud  (within  several  hours  of 
detonation)  particles  represent  both  nx  and  n2  behavior  as  indicated  by  Nathans’ 
analysis  of  surface  burst  clouds  (16,64).  As  should  be  expected,  close-in  ground  sam¬ 
ple  particles  are  more  characteristic  of  n2  as  noted  by  Tompkins  (predominantly 
irregular  in  shape,  ref.  49:396)  and  Norment  (DELFIC  distribution  more  closely 
approximates  n2,  ref.  8).  Air  bursts  with  little  or  no  surface  material  present  such 
that  all  debris  reaches  fireball  temperatures  produce  an  rjj  type  distribution  almost 
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without  exception  (exceptions  include  some  tower  shots,  and  lower  air  bursts  with 
some  limited  surface  material  contamination).  A  heuristic  illustration  of  the  origin  of 
the  nx  and  n2  distributions  appear*  in  figure  53. 

The  n1  parameters  are  well  established  by  the  independent  results  previously 
discurjed-  namely,  a  lognormal  distribution  with  rm  of  .1  n  and  slope  of  2  describes 
the  behavior  of  stratospheric  tracers,  air  burst  debris,  and  unmixed  fission  products. 
The  establishment  and  corroboration  of  n2  parameters  was  not  as  straightforward. 

£.  Characterization  of  n2. 

The  most  direct  evidence  available  for  an  actual  n2  size  distribution  is  Nathans’ 
analysis  of  cloud  samples  from  the  10  KT  silicate  surface  burst  discussed  in  section 
V.C.2.  Nathans  divided  his  samples  into  spherical  and  irregular  particle  fractions 
and  characterized  their  size  distributions  separately.  The  results  for  one  sample  are 
shown  in  figure  51.  His  n2  distribution  is  lognormal  with  a  median  radius  of  .17  /i 
and  a  slope  of  2.  Unfortunately,  this  is  the  only  instance  in  which  he  intentionally 
characterized  spherical  and  irregular  particle  size  distributions  separately  (68). 

Nathans  did  characterize  what  is  purported  to  be  the  undifferentiated  size  distri¬ 
butions  of  clouds  from  four  surface  bursts  (Johnie  Boy,  Koon,  Zuni,  and  Bravo,  ref. 
18).  However,  a  careful  reading  of  his  documentation  reveals  that  (for  undisclosed 
reasons)  he  included  only  irregular  particles  in  his  size  analysis  (16:367).  Thus,  the 
particle  size  statistics  for  these  events  were  used  to  get  additional  information  on  n2. 

A  puzzling  aspect  of  Nathans’  particle  statistics  for  the  four  surface  bursts  was 
the  power  law  behavior  (not  lognormal)  of  the  size  distributions  above  1/i  diameter 
(see  figures  4-6).  This  behavior  can  be  explained  by  noting  the  position  of  the  air¬ 
craft  with  respect  to  the  cloud  at  the  time  of  sampling.  Apparently  the  sampling  air¬ 
craft  could  not  get  up  to  the  cloud  altitude.  In  all  but  two  of  Nathans’  samples  for 
the  four  events,  samples  were  extracted  at  altitudes  below  the  visible  cloud  (see  table 
V).  Only  Johnie  Boy  samples  827L  and  842L,  although  lower  than  the  cloud  vertical 
centroid,  were  within  the  visible  cloud.  It  was  apparent  that  the  low  sampling  alti¬ 
tudes  affected  the  size  distribution  of  the  samples  vis-a-vis  the  aggregate  cloud. 

Nathans  did  correct  for  sedimentation  in  plotting  his  size  distributions  (16:379). 
However,  his  correction  factor  did  not  account  for  a  variable  vertical  cloud 
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4  hours  post  shot 


Figure  54.  Group  Fall  Illustration 
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Table  V 

Nathans’  Samples:  Bravo,  Zuni,  Koon,  Johnie  Boy 

Event/sample 

Yield 

Soil 

Cloud 

Top*  Bottom 

Sample 

Time  Altitude 

Bravo  (1954) 

15  MT 

coral 

34km 

16.6km 

4hrs 

15.5km 

Zuni  (1956) 

3.4  MT 

coral 

24.5km 

15km 

S.lhrs 

12.7km 

Koon  (1954) 

150  KT 

coral 

16.5km 

/1086 

3.1hrs 

11.95km 

/051 

3.1hrs 

11.95km 

/7269 

3.1hrs 

11.95km 

Johnie  Boy  (1962) 

.5KT 

silicate 

5.18km 

3.81km 

/842R 

.6hrs 

2.95km 

/842L 

.4hre 

3.4km 

/245L 

.8hre 

3.7km 

/827L 

^^5hrs^ 

4.35km 

*In  fallout  calculations,  the  cloud  centroid  was  taken  as  the  mean  of  top  and  bottom 
altitudes.  The  stabilized  cloud  was  modeled  as  Gaussian  in  the  vertical  dimension  with 
at  -  \/4(CT-0B) 

*A  strong  wind  shear  at  about  11km  directed  the  bottom  portion  of  the  cloud  away 
from  the  remainder  of  the  cloud. 


concentration  profile  and  consequently  was  not  very  large  (less  than  a  factor  of  2  in 
all  cases).  Thus,  his  corrected  distributions  are  very  nearly  equal  to  what  was  actu¬ 
ally  measured.  Nathans’  correction  factor  does  account  for  the  vertical  velocity  gra¬ 
dient  which  does  not  produce  as  large  an  effect  as  the  vertical  concentration  gradient 


1  (Section  A.3.c). 

Nathans  used  aircraft  samples  taken  between  .4-4  hours  after  detonation.  Fig¬ 
ure  51  illustrates  size  group  stratification  predicted  (assuming  Gaussian  dispersion 
and  using  subroutine  STRATFAL)  at  4  hours  following  the  Bravo  event.  The 
graph  is  plotted  with  n(r)  as  constant  so  that  only  variations  due  to  the  vertical  con¬ 
centration  gradient  are  depicted.  Note  that  at  any  sampling  altitude,  ze  (a  12  KM 
sampling  altitude  is  marked),  the  sample  n(r)  is  greatly  affected  by  the  vertical 
stratification.  At  time  t=0,  assuming  that  all  size  groups  are  distributed  uniformly 
in  the  z  direction,  a  distribution,  n0(rj ,t  asO,zi ),  proportional  to  the  true  cloud  distri¬ 
bution  would  be  measured  at  any  altitude.  If  a  sample  is  taken  at  a  later  time,  a  dis¬ 
tribution  of  nt(r},tt ,zt )  is  measured.  It  is  then  possible  to  relate  ng  to  n,  by  using 
the  correction  factor,  Cj\ 
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Cj  is  defined  as 


n0{rrO,  zt) 


t,  za) 


(63) 
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where  f  is  the  vertical  distribution  function.  In  this  case,  analogous  to  equation  30, 
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f(z  -  z»)  - 


-% 


V27T  <7t 


i-i. 


(65) 


Cj  fur  Bravo  is  plotted  in  figure  55  (DELF1C  fall  mechanics  were  used).  Note  that 
the  correction  factor  peaks  at  that  particle  diameter  which  has  just  reached  aircraft 
altitude  at  sampling  time.  The  apparent  population  at  this  diameter  would  be 
enhanced  in  the  sample.  The  sample  would  be  depleted  in  larger  particles  since  they 
would  have  fallen  past  the  sampling  altitude.  The  smallest  particles  would  be  present 
in  their  true  proportions  (Cj  =  1). 
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Figure  56.  Size  Distribution  Correction  Factor 
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It  is  suspected  that  this  correction  factor  boosted  the  population  in  Nathans’ 
samples  in  the  10-  100m  range  and  that  this  effect  causes  a  lognormal  distribution  to 
approximate  a  power  law.  This  possibility  was  investigated  by  using  Nathans’  log¬ 
normal  characterization  of  irregular  particles  for  the  10  KT  silicate  surface  burst  as 
the  aggregate  cloud  n2  distribution  for  the  four  events.  For  the  10  KT  event, 
Nathans  realized  that  sedimentation  had  affected  the  population  of  particles  larger 
than  lO/i  in  his  samples  and  therefore  confined  his  analysis  to  particles  of  diameter  2/x 
and  smaller  (see  figure  52).  Gravity  sorting  would  not  have  significantly  affected  the 
relative  population  of  these  small  particles  at  sampling  time  (figure  10).  Thus  the  10 
KT  particle  statistics  were  a  reasonable  standard  of  comparison.  It  was  suspected 
that  the  lognormal  parameters  developed  for  <2/i  particles  would  also  describe  the 
larger  particle  population. 

Starting  with  a  lognormal  aggregate  cloud  distribution,  the  correction  factor  was 
applied  to  determine  the  sample  population  at  the  time  and  altitude  of  extraction. 
The  calculation  for  Zuni  sample  049  is  illustrated  in  figure  56.  Application  of  the 
correction  factor  to  the  lognormal  distribution  yields  a  curve  that  does  in  fact 
approximate  a  power  law  function.  Indeed,  a  straight  line  fit  to  the  curve  falls  as 
r-39  (Nathans’  fit  yielded  an  exponent  of  -4.07).  Similar  results  are  obtained  for  the 
other  samples.  In  each  case,  the  correction  factor  applied  to  a  lognormal  distribution 
with  rm  ass  .17  /x  and  0  =  ln(3)  masquerades  as  a  power  law  function.  In  order  to 
determine  the  best  slope,  rm  was  held  constant  and  0  was  varied  until  the  slope  of 
the  corrected  lognormal  distribution  matched  Nathans’  power  iaw  exponent.  In  each 
case  a  logarithmic  slope  of  '3  provided  the  best  fit  (see  figures  57-59).  Varying  rm 
between  .1— *.5/i  did  not  significantly  affect  the  optimum  0. 

The  correction  factor  as  expressed  in  equation  64  did  not  explain  the  depletion  of 
submicron  particles  in  Nathans’  samples  (figures  4-6).  A  more  elaborate  correction, 
taking  into  account  variations  in  f(z)  with  size  group  is  probably  needed.  Since 
Nathans  did  not  find  spherical  particles  in  samples  taken  underneath  the  visible 
cloud,  there  is  reason  to  believe  that  smaller  particles  were  more  tightly  clustered  at 
the  cloud  center  (smaller  <7t)  such  that  their  presence  at  sampling  altitude  would  have 
been  considerably  diminished  (65:3).  Such  clustering  could  be  caused  by  the  centrifu¬ 
gal  action  of  the  toroidal  cloud  motion  which  should  preferentially  disperse  the 


Figure  57.  Lognormal  vs.  Power  Law  Fit, 
Johnie  Boy  Sample  245L 


Figure  58.  Lognormal  vs.  Power  Law  Fit, 
Johnle  Boy  Sample  842L 


Figure  59.  Lognormal  vs.  Power  Law 
Fit,  Bravo  Sample 
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heavier  particles.  This  phenomenon  would  explain  the  absence  of  spherical  particles 
in  Bravo,  Koon,  and  Zuni  samples  (all  significantly  below  the  cloud  centroid)  and  the 
presence  of  spherical  particles  in  Johnie  Boy  samples  (taken  near  the  could  centroid). 
Another  possible  explanation  for  small  particle  depletion  is  agglomeration  due  to 
differences  in  group  fall  velocity.  The  relatively  motionless  small  particles  might  have 
been  depleted  by  sticking  to  larger  particles  falling  through  their  midst. 

Admittedly  there  are  many  remaining  questions  concerning  the  n2  distribution. 
Only  one  instance  (10  KT  silicate  burst)  was  found  in  which  n2  was  intentionally 
analyzed  as  a  separate,  independent  entity.  It  appears  that  only  n2  type  particles 
were  analyzed  for  Johnie  Boy,  Bravo,  Zuni  and  Koon.  It  is  puzzling  that  few  if  any 
spherical  particles  were  present  in  the  Bravo,  Zuni,  and  Koon  samples.  A  possible 
explanation  has  been  offered  but  there  may  be  others.  For  instance,  Pacific  events 
may  not  have  produced  significant  numbers  of  spherical  particles  (although  some  were 
observed  on  the  ground  near  Tewa).  Even  though  present  results  indicate  a  fairly 
consistent  n2  parameter  values  among  events,  significant  excursions  are  expected  since 
these  particles  originate  from  surface  materials  which  differ  considerably  in  morphol¬ 
ogy  and  thermal  history  (nx  originates  from  fireball  vapors  whose  initial  temperature 
is  nearly  independent  of  surface  material  and  yield).  Indeed  there  is  evidence  that 
for  buried  bursts  n2  is  multimodal  (60).  And  Norment’s  nominal  DELFIC  distribu¬ 
tion  implies  that  larger  logarithmic  slopes  exist. 

Nathans’  10  KT  silicate  surface  burst  results  provide  the  most  solid  evidence  for 
the  existence  of  an  independent  n2  distribution.  The  10  KT  result?  are  consistent 
with  Nathans’  distributions  for  Johnie  Boy,  Bravo,  Zuni,  and  Koon  if  a  Gaussian 
cloud  stratification  correction  factor  is  applied.  Assuming  these  n2  parameters  to  be 
representative  of  clouds,  an  important  remaining  question  is  what  are  the  relative 
populations  of  rij  and  n2?  The  relative  number  (or  mass)  of  particles  in  each  distri¬ 
bution  determines  the  optical  properties  and  longevity  of  nuclear  clouds. 

F.  The  nj  —  n2  Split. 

The  measured  r»j/n2  ratio  for  cloud  samples  from  the  two  silicate  surface  bursts 
discussed  in  section  V.E  were  nearly  identical.  The  ratio  was  2.3  for  Johnie  Boy  and 
2.2  for  the  10  KT  silicate  surface  burst.  These  ratios  would  not  have  been 
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significantly  affected  by  sedimentation  since  nearly  100%  of  the  population  lies  below 
5 fi  for  both  and  n2.  The  mass  ratio  may  then  determined  from: 
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(86) 


assuming  particle  densities  (plt  p2)  are  roughly  the  same  between  the  modes.  Using  a 
lognormal  n,  with  rm  of  .1  n  and  slope  of  2,  and  a  lognormal  n2  with  rm  of  .17  p 
and  slope  of  3,  and  Nx/N2  *»  2.2,  the  result  is: 
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(67) 


Thus,  even  though  n  i  has  a  higher  population,  it  represents  only  1.7%  of  the  mass 
lofted.  The  mass  ratio  expressed  in  equation  67  can  be  used  to  compute  the  com¬ 
bined  fraction  of  mass  below  lp,  /m,  which  equals  5.8%.  This  result  is  close  to  the 
TTAPS  value  of  8.4%. 


Since  the  mass  of  soil  lofted  decreases  with  height  of  burst  (13:3-7),  Mx/M2 
should  increase  accordingly.  An  indication  of  increase  with  HOB  may  be  obtained  by 
estimating  the  fraction  of  activity  carried  by  nj  vs  n2.  Local  fallout  data  can  be  used 
to  get  an  indication  of  the  activity  split.  Davis  (54)  has  attempted  to  derive  size  dis¬ 
tributions  from  local  activity  deposited  vs  time.  Unable  to  model  air  burst  deposition 
using  a  unimodal  lognormal  distribution,  he  found  that  a  bimodal  distribution  couid 
explain  the  data.  He  fit  the  observed  rate  of  local  activity  deposition  for  actual 
events  of  varying  scaled  height  of  burst  using  the  following  expression: 


oo  oo 

Fx  J  Oj(r )dr  +  F2  J  a2{r)dr 

r(t) 


(68) 


where  A0  is  the  total  activity  lofted,  r(t)  is  the  radius  of  particles  grounded  at  time  t, 
F 2  is  the  fraction  of  activity  carried  by  n2,  and  Fx  —  1  —  F 2.  By  fixing  nj  (lognor¬ 
mal  with  .  1/i  and  slope=2)  and  assigning  n2  an  rm  of  .2  /j  Davis  obtained  the 
results  in  Table  VI  (plotted  in  figure  60).  The  results  indicate  that  10  -  100%  of  the 
activity  is  carried  by  n2  for  a  <3  foot  scaled  height  of  burst  (SHOB=HOB*/CT’_1/3) 
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Figure  60.  Fraction  of  Activity  on  i\2  vs.  Scaled  Height  of  Burst 


but  that  the  fraction  drops  rapidly  above  100  scaled  feet.  If  the  activity  is  distri¬ 
buted  as  the  mass  of  particles,  then  51(r)  and  S2(r)  are  constant  and  the  mass  ratio 
is  proportional  to  the  activity  ratio: 


£1 


OO 

4-*Pi/ni(»VS|(r)dr 

*  0 _ 

OO 


-i*rp2/n2(r)r3S2(r)dr 


<Si>  Mx 

<52>  M2 


(69) 


Davis’  Fx/F 2  vs  scaled  height  of  burst  is  plotted  in  figure  61.  These  results  will  of 
course  change  with  assumed  particle  size  parameters  and  should  be  used  with  caution. 
However,  the  indicated  trends  are  correct.  It  is  interesting  that  the  functional  form 
of  F 2  vs  SHOB  (figure  80)  is  reminiscent  of  Carpenter’s  fit  to  the  total  mass  lofted  vs 
SHOB  (ref.  70,  figure  62).  This  is  logical  since,  for  higher  bursts,  less  surface  mass  is 
available  to  mix  with  fission  products,  hence  f2  should  decrease  roughly  proportional 
to  mass  lofted.  It  is  interesting  that  Davis’  cliff  in  mass  lofted  occurs  at  ~200  scaled 
feet,  noticeably  lower  than  the  '700  ft.  cliff  exhibited  by  Carpenter’s  fit. 

Notice  that  if  the  activity  is  radially  distribmed  on  particles  it  is  not  necessary 
to  know  specific  activities,  but  rather  mass  ratios  can  be  determined  from  number 
ratios  via  simple  scaling  relations: 

Ml  -di/<^i>  iVj  <rp> 

M2  "  A2/<52>  "  JV2  <r3>  '  0) 


These  number  ratios  can  be  readily  determined  from  microscopic  examination  of  fal¬ 
lout  samples  preserved  from  atmospheric  tests.  If  3ize  distributions  are  known  then 
average  specific  activity  ratios  may  be  determined  from  a  knowledge  of  the  total 
activity  carried  by  each  mode.  Aggregate  area  fractions  also  obey  a  simple  scaling 
relation: 


Area!  M,  <r2>  <r,2> 

- L  wm — 1 - : - 1 -  (71) 

Areo2  M2  <r,3>  <r22  > 

The  lognormal  behavior  of  r»  2  and  n2  further  simplifies  the  moments  calculations. 

The  results  in  Table  VI  give  evidence  of  the  suspected  variability  of  n2.  The 
logarithmic  slope  varies  from  3.2  to  6,  with  a  mean  of  4.2  wh'ch  is  near  the  DELFIC 
nominal  value  of  4.0.  While  the  uncertainties  inherent  in  Davis’  calculations  are  large 
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Fiqure  62.  Oust  Mass  Lofted  into  the  Stabilized  Cloud  (Showing  Estimated  Fraction  of  Glass 
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Table  VI 


Lognormal  Parameters:  Bimodal  Distribution8 


Event 

F 2  Median  Radius  rm# 

Log  Slope  (e^*) 

Surface  Bursts 

PLUMBBOB  Coulomb-B 

0.214 

0.2 

4.7 

BUSTER-JANGLE  Sugar 

0.924 

0.2 

4.6 

Tower  Bursts 

TEAPOT  Hornet 

0.115 

0.2 

5.2 

TEAPOT  Bee  . 

0.0282 

0.2 

6.0 

TEAPOT  Apple 

0.0404 

0.2 

4.4 

PLUMBBOB  Shasta 

0.0994 

0.2 

4.8 

UPSHOT-KNOTHOLE  Nancy 

0.0871 

0.2 

4.3 

TEAPOT  Apple  II 

0.0852 

0.2 

4.1 

UPSHOT-KNOTHOLE  Harry 

0.703 

0.2 

3.4 

UPSHOT-KNOTHOLE  Simon 

0.195 

0.2 

3.6 

Air  Bursts 

PLUMBBOB  Morgan 

0.00583 

0.2 

4.7 

PLUMBBOB  Owens 

0.0407 

0.2 

3.2 

UPSHOT-KNOTHOLE  Grable 

0.0240 

0.2 

3.4 

PLUMBBOB  Stokes 

0.00036 

0.2 

4.1 

PLUMBBOB  Priscilla 

0.0366 

0.2 

3.3 

PLUMBBOB  Hood 

0.00212 

0.2 

4.0 

4  Nominal  nt  assumed  with  rm|  —  .l/i  and  (3X  —  ln2 
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and  difficult  to  quantify,  his  results  give  further  support  to  the  existence  of  a  bimodal 
particle  distribution  in  clouds  and  also  give  a  rough  indication  of  activity  partitioning 
between  the  modes. 
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VI.  Summary  and  Conclusions. 


This  chapter  will  discuss  the  ramifications  of  Chapter  V  results  and  suggest 
further  applications  of  atmospheric  test  data  to  the  nuclear  winter  problem. 


A.  Blmodal  Cloud  Debris. 

The  m^jor  conclusion  of  this  study  is  that  cloud  debris  behavior  can  be 
explained  by  the  existence  of  a  bimodal  size  distribution.  The  puzzling  dichotomy 
between  early  time  and  late  time  debris  sedimentation  (vividly  demonstrated  by  the 
failure  of  the  stratospheric  size  distribution  to  produce  appreciable  intermediate/ local 
fallout  by  the  sedimentation  process)  can  be  explained  if: 


n(r)  -  n^r)  +  n2(r) 


I  *  > 

V2 *0, 


— Vfc 


N2e 


ln[r)  -  in(rmt)  f 

A  I 


V27T  j32 


(72) 


where  n  j  governs  the  late  time  behavior  of  activity  removal  and  n2  governs  the  early 
time  behavior.  Figure  63  plots  n(r)  as  expressed  above  using  Nathans’  measured  10 
KT  silicate  surface  burst  parameters.  The  TTAPS  and  DELFIO  distributions  are 
included  for  comparison.  Supporting  evidence  for  the  bimodai  distribution  includes 
the  size  distribution  indicated  by  the  bulk  vertical  transfer  of  stratospheric  tracers 
over  a  period  of  months  to  years  (section  VA)  and  Nathans’  analysis  of  early  time 
cloud  debris  from  surface  bursts  (section  V.C).  The  n  t  (or  “hot'*)  component  of  cloud 
debris  appears  to  originate  from  the  condensation  of  material  vaporized  in  the  fireball 
and  is  present  for  both  surface  and  air  bursts.  The  n2  (or  "cold")  component  appears 
to  be  unmelted  or  partially  melted  surface  material  which  never  reached  fireball  tem¬ 
peratures.  The  n2  component  is  not  present  in  cloud  samples  from  free  air  bursts. 
Table  VII  contrasts  the  properties  of  and  n2. 

The  microscopically  measured  size  distribution  parameters  for  the  n  j  component 
are  fairly  consistent  from  event  to  event  (17).  It  is  significant  that  the  parameters 
properly  predicting  the  removal  of  activity  from  the  stratosphere  are  consistent  with 
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Table  VH 


Characteristics  of  Nx\a  N2  Particles 


Nx 

spherical 

irregular 

high  specific  activity 

low  specific  activity 

5,(r)  a  r3 

S^r)  a  r2 

Origin: 

Origin: 

vaporized  material 

unmelted/partially  melted  material 

weapon  mass  -+-  soil 

soil 

Size  Distribution: 

Size  Distribution: 

lognormal 

lognormal 

rm  » 

rm  ~  -2 fi 

0t*ln{  2) 

0  ln( 3)  —  /n(5) 

Stratospheric  Removal: 

Stratospheric  Removai: 

Sedimentation  +  Diffusion 

Sedimentation 

PRESENCE*: 

|  Air  burst . yes 

. no 

Low  air/tower... varies  with  HOB 

. varies  with  HOB 

Surface... small  fraction  of  mass 

...large  fraction  of  mass 

‘Reference  76 


RELATIVE  FREQUENCY 


Figure  63.  Bimodal  Distribution  Compared  to  Nathans'  and  DELFIC 


the  microscopically  measured  size  distributions  for  hot  particles  from  air  and  surface 
bursts  (figures  35,30,42).  Representative  values  for  the  median  radiu3  and  slope  are 
.1/i  and  2  respectively.  These  values  predict  the  average  behavior  of  stratospheric 
**Sr  debris  injected  by  97  events  over  a  15  year  period. 

Characterization  of  n2  has  been  less  satisfactory  due  to  data  limitations  and  the 
expected  variation  in  n2  parameters  with  surface  material  and  burst  height.  Analysis 
of  limited  data  for  surface  bursts  provided  by  Nathans  indicates  a  median  radius  and 
slope  in  the  vicinity  of  .17/u  and  3  for  the  n2  component.  Davis’  results  (54)  indicate 
that  logarithmic  slopes  may  vary  considerably  about  a  mean  of  ~4.  Clearly,  more 
work  'is  needed  to  bound  the  size  and  mass  properties  of  the  "cold"  distribution.  It 
may  be  that  n2  is  itself  multimodal,  especially  for  buried  bursts  (69).  In  addition 
the  relation  of  n2  fraction  to  burst  height  needs  to  be  resolved  with  confidence.  This 
relation  is  bound  up  with  the  problem  of  determining  mass  lofted  per  kiloton. 


The  existence  of  a  bimodal  cloud  distribution  leads  to  a  new  equation  for  com¬ 
puting  the  total  mass  lofted: 


Mf  *■  Af 


<sx>  <s2> 


(73) 


The  present  enorc  not  investigated  the  full  implications  of  this  equation  which 
enables  computation  of  mass  lofted  based  upon  a  knowledge  of  total  activity  depo¬ 
sited.  Reasonable  estimates  of  Fx  and  F2  should  be  possible  from  a  careful  evalua¬ 
tion  of  existing  fallout  data.  For  example,  Tompkins  (71)  has  developed  fits  for  the 
fraction  of  activity  deposited  locally  based  upon  data  from  several  Nevada  events. 

|  Since  a  large  part  of  the  n2  distribution  and  almost  none  of  the  n  j  distribution  falls 

|  in  24  hours,  Tompkins’  fraction  roughly  approximates  F2,  at  least  for  above  ground 

I 

shots. 

There  is  considerable  uncertainty  in  the  average  specific  activity  parameters. 
The  data  required  to  narrow  these  uncertainties  does  not  appear  to  be  plentiful,  par- 
|  ticularly  for  determining  the  variation  of  <5j>  and  <S2>  with  height  of 

burst.  However  data  from  surface  bursts  and  free  air  bursts  ure  available,  which 
taken  together  with  some  reasonable  physics  could  be  used  to  develop  a  credible 
model.  For  purposes  of  illustration,  consider  an  approximate  calculation  for  Johnie 
I  Boy  whose  fission  yield  was  .5  KT  (equivalent  to  ~7  x  1022  fissions).  Tompkins  (71) 


^ 
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indicates  that  "70%  of  the  activity  from  Johnie  Boy  was  deposited  locally.  Using 
this  as  a  rough  value  for  and  estimating  <St>  and  <S2>  to  be  ~10!6  and  "1015 
respectively,  then  MT  «  .051  AT.  This  translates  to  .12  kilotons  of  dust  lofted  per 
kiloton  of  yield. 

B.  Cloud  Stratification  Effects. 

An  important  conclusion  from  the  investigation  of  the  power  tail  behavior  of 
sample  from  Johnie  Boy,  Bravo,  Zuni  and  Koon  (section  V.E)  is  that  size  distribu¬ 
tions  vary  tremendously  with  sample  location  and  time  due  to  sedimentation  induced 
vertical  cloud  stratification.  Consequently,  extreme  caution  should  be  exercised  in 
generalizing  a  cloud  distribution  from  the  distributions  of  isolated  samples.  The 
results  show  specifically  that  clouds  whose  total  particle  ensemble  is  lognormally  dis¬ 
tributed  can  yield  sample  distributions  which  obey  a  power  law  function.  It  is  quite 
possible  that  Natl  .ns’  power  law  distributions  originated  from  clouds  whose  aggre¬ 
gate  population  was  closer  to  lognormal. 

Cloud  stratification  effects  on  sample  distributions  become  extremely  pronounced 
with  increasing  distance  from  the  cloud  centroid.  The  stratification  effect  may  in  fact 
explain  the  narrow  distributions  observed  for  ground  samples  which  seem  to  vary  so 
much  with  distance  downwind.  An  example  of  the  size  distribution  which  would  be 
measured  on  the  ground  12  hours  post  shot  underneath  a  Gaussian  cloud  with  a  stan¬ 
dard  devia.  on  of  700m  is  shown  in  figure  64.  The  distribution  peaks  at  35 n  with  a 
half-width  of  only  5u.  The  half-width  is  a  function  of  the  sedimentation  correction 
factor  (equation  61)  which  »s  in  turn  a  function  of  the  cloud  vertical  density  distribu¬ 
tion.  Thus,  it  may  be  possible  to  learn  something  about  the  vertical  density  of  clouds 
by  studying  the  half  width  of  particle  size  distributions  of  ground  samples. 


C.  Optical  Depth  Implications. 

Using  the  bimodal  size  distribution  expressed  in  equation  73,  it  is  possible  to 
compute  the  optical  thickness  vs  time  for  a  stratospheric  dust  cloud  representative  of 
the  TTAPS  5000  Megaton  exchange.  Analogous  to  equation  10: 


OT 


3Mp  Qe 


<r2> 


<r3> 


(74) 
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Figure  64.  Ground  Sample  Olstributi 
12  Hours  Post-Shot 


Using  equations  6  and  7  for  the  2nd  and  3rd  moments  of  the  lognormal  iistributions 
n  j  and  n2  the  following  expression  for  the  initial  optical  thickness  results 


OT 


3  MtQ, 
4ptAe 


.  ,2/n(r„+2^D 
JV2 


3/n(r.,+l^D 
+  e  2 


(75) 


Using  the  cloud  parameter  values  from  section  I.B.1 
(Ae  m  1.14  x  1014  m2,  pt  —  2600  kg/m3,  and  Qt  —  2.5)  the  initial  optical  thickness 
is  2.0.  This  is  comparable  to  the  initial  dust  optical  thickness  of  1.25  computed  by 
the  TTAPS  group.  Recall  (figure  3)  that  the  nominal  DELFIC  n(r)  predicts  an  initial 
O.T.  of  ,25,  and  the  Ramaswamy-Kiehl  (R-K)  distribution  predicts  11.6.  The 
bi modal  distribution  was  used  in  subroutine  STRATFAL  to  predict  the  rate  of  decay 
of  optical  thickness  for  a  cloud  initially  centered  at  25  KM.  The  results  are  overlaid 
with  DELFIC,  TTAPS,  and  R-K  in  figure  65.  The  bimodal  distribution  optical 
attenuation  is  initially  higher  than  the  TTAPS  result,  but  decays  more  rapidly  level¬ 
ing  off  to  a  value  of  of  the  TTAPS  optical  thickness  (the  factor  of  Vi  is  certainly 
within  model  uncertainties).  The  curves  are  parallel  at  late  times  because  the  remain¬ 
ing  particles  are  small  enough  that  their  removal  is  governed  by  the  exponential  tur¬ 
bulent  diffusion  effect.  Thus,  the  size  distribution  derived  in  this  research,  although 
different  than  the  size  distribution  used  in  the  TTAPS  study,  yields  similar  estimates 
of  optica!  thickness  magnitude  and  decay  (other  effects  being  equal).  Although  the 
power  law  behavior  may  well  be  an  aberration  of  cloud  stratification,  Nathans’  distri¬ 
bution  appears  to  provide  a  reasonable  estimate  of  sunlight  attenuation  magnitude 
and  duration.  The  similarity  of  the  results  is  due  to  the  fact  that  the  submicron  frac¬ 
tions  are  comparable. 


D.  Recommendations. 

1.  The  evidence  supports  the  existence  of  a  bimodal  particle  size  distribution  in 
nuclear  clouds.  The  distribution  expressed  in  equation  72  should  be  seriously  con¬ 
sidered  for  predicting  optical  attenuation  effects  and  fallout  removal.  The  existence 
of  distributions  of  different  consistencies  and  thermal  origins  (Table  VII)  warrants  an 
investigation  of  possibly  large  differences  in  optical  properties  as  well. 
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Figure  65.  OT  vs.  Time,  Cloud  Initially  Centered  at  25  km 
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2.  Work  is  needed  to  characterize  n2(r).  This  research  has  uncovered  only  one 
instance  of  separate  size  analysis  of  spherical  and  irregular  particles.  There  may  be 
others.  If  not,  it  is  recommended  that  further  analysis  of  samples  from  atmospheric 
test  clouds  be  conducted  and  that  any  such  analysis  develop  separate  statistics  for 
spherical  and  irregular  particles.  The  n2  parameter  values  used  here  (median  radius 
A7fi  and  slope  of  3)  are  based  on  analysis  of  actual  data  by  Marcel  Nathans,  and  are 
certainly  reasonable  values  to  use  in  the  interim  for  estimates  of  surface  burst  effects. 
Indications  are  that  the  slope  ranges  higher.  Thus  the  nominal  DELFIC  distribution 
is  also  a  plausible  representation  of  n2. 

3.  Computations  of  mass  lofted  per  kiloton  should  be  reexamined  in  light  of 
the  bimodal  nature  of  the  cloud  distribution.  An  attempt  should  be  made  to  quan¬ 
tify  the  parameters  of  equation  73  based  upon  test  data  over  a  range  of  scaled  burst 
height/depth.  It  is  reasonable  to  expect  that  the  uncertainty  in  estimates  of  mass 
lofted  vs  yield  and  height  of  burst  can  be  considerably  reduced. 

4.  It  would  be  useful  to  run  the  rigorous  models  of  climate  effects  from  nuclear 
war  using  the  bimodal  n(r)  function  and  parameters  developed  in  the  present 
research.  Although  it  appears  that  optical  effects  are  not  drastically  different  from 
TTAPS,  the  effect  on  temperature  drop  needs  to  be  investigated.  The  higher  initial 
magnitude  and  more  rapid  decay  rate  may  cause  larger  differences  in  the  predicted 
temperature  changes. 
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APPENDIX  A 


TREATMENT  OF  FALL  MECHANICS 


Fall  mechanics  involves  the  solution  of  the  one  dimensional  force  balance  equa¬ 
tion  including  a  viscosity  term: 

(A.1) 

where  m  is  the  particle's  mass,  u  is  its  velocity,  and  /  is  the  viscous  damping 
coefficient.  Particles  are  assumed  to  be  spherical.  The  functional  relationship 
between  the  viscous  force,  ftt ,  and  particle  diameter  varies  with  particle  size  and  alti¬ 
tude.  At  sea  level,  for  particles  less  than  10/i,  Stoke’s  law  applies: 


/  —  8  irqr 


(A.2) 


and  the  force  balance  equation  becomes: 


6  Tr\r  —  -2-tt  rzpg 
o 


(A.3) 


where  77  is  the  dynamic  viscosity  (kg/sec/m)  and  r  is  the  particle’s  radius.  Note  that 
the  resisting  force  is  proportional  to  radius.  Solving  for  terminal  velocity: 

2  r2pg 


uterm  ™ 


Qtj 


(A.4) 


For  particles  greater  than  aerodynamic  drag  is  significant  and  the  force 

balance  equation  becomes: 

Vtpt  u 2  Cd  7r r2  -  jnr3pg  (A.5) 


The  solution  of  this  equation  is  complicated  by  the  fact  that  Cd  is  a  function  of  the 
particle’s  velocity.  McDonald  (73)  suggested  a  method  of  solution  by  defining  u  in 
terms  of  the  Reynolds  number: 


u  » 


Rev 

%Par 


(A.8) 
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McDonald  then  defined  a  quantity,  Q : 


32 pa(p  -pa)gf3 
3  rf- 


(A.7) 


Q 


R?cd  - 


which  Davies  (74)  has  fit  to  various  polynomial  expressions  depending  upon  its  mag¬ 
nitude.  Q  is  referred  to  as  the'Davies  number").  For  a  given  altitude,  the  fall  velo¬ 
city  is  computed  by  calculating  Q ,  finding  Rt  from  Davies’  polynomial  fits,  and  then 
computing  it  from  Re .  Polynomial  fits  for  determining  Re  are  as  follows: 


for  140  <  Q  <  4.5E7: 

log l0Re  -  -1.29538  +  986log10Q  +  .O466771og20<3  +  1.1235£  -Slog&G  (A.8) 


for  84.175  <  Q  <  140: 

Re  -  Q(  4.166667E-2  -  2.3383£-4^  +  2.0154£-8<?2  -  6.9105<?3>  (A.9) 


for  .3261  <  Q  <  84.175: 

Re  -exp[-3.18657+.992696/nG-.153193£-2/n2<?-.987059£~3/n3<3-  578876£-3/n4g 

+.855l76£-4/n5g-.327815£-5/n6Qj  (A.10) 


for  Q  <  .3261 


Stokes  equation  (A.5)  applies. 

For  particle  diameters  approaching  the  mean  free  path  (L)  of  air  molecules,  the 
drag  for  a  given  velocity  is  less  than  predicted  by  the  polynomial  fits  and  a  slip 
correction  factor  must  be  introduced.  For  example,  when  d/L=10  the  Stokes  equa¬ 
tion  underestimates  u  by  about  15  percent  (9).  At  50,000  ft,  L=.4  microns,  so  it  is 
clear  that  the  slip  factor  is  important  for  gravity  fall  calculations. 


The  Cunningham  slip  correction  factor  used  in  HASP  analysis  (10)  and  DELFIC 
(8)  was  derived  for  d/L  ratios  between  2  and  100.  Knudsen  and  Weber  (57)  proposed 
a  general  formula  that  fits  data  over  a  wider  range  of  d/L  values: 


C 


+  A  2exp 


-2A3r 


(A.11) 


130 


■where  A^l.644,  A2=.552,  and  A3=.656.  The  Knudsen-  Weber  slip  correction  for¬ 
mula  was  incorporated  into  the  fall  velocity  model. 

The  equations  in  this  appendix  were  encoded  as  shown  in  the  listing  which  fol¬ 
lows.  Computed  terminal  velocities  compare  to  within  20%  of  those  computed  by 
models  used  in  the  High  Altitude  Sampling  Program  for  particle  diameters  ranging 
from  .001  -10/*  (10:154). 
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FUNCTION  VEL ( ZB , RP ) 

C  COMPUTES  TERMINAL  VELOCITY  OF  SPHERICAL  PARTICLES  USING 

-FULL  DELFIC 

FALL  MECHANICS  AND  U.S.  STANDARD  ATMOSPHERE  EQUATIONS 
SUBROUTINE  ATMOS  RETURNS  DENSITY ,  VISCOSITY  AND  MEAN  FREE 
PATH  AT  ALT  ZB 

REAL  MFP , LOORBY 
C  EXTERNAL  ATMOS 

CALL  ATMOS(ZB,DENS, VISC, MFP) 

DENSF«2600. 

0*9.8 

QZ*4*DENS*(DBNSF-DENS)*G*(2 .*RP)**3/(3.*VISC*VISC) 

C  WRITE(6, 1051 )QZ, DENS, VISC, MFP 

Cl  051  FORMAT ( 'QZ, DENS, VISC, MFP' , 4 ( B 1 0 . 3 ) ) 

YYsALOG(QZ) 

C  DENSF  IS  PARTICLE  DENSITY,  Q  IS  GRAV ,  QZ  IS  DAVIES  NUMBER 

DEFINED  IN 

C  DNA  TR-5159F-1 ,  P24 

IP ( QZ  .LE.  . 326 1 ) THEN 

VEL»VISC*QZ/(24.*DEN3*2.*RP) 

GO  TO  1100 
ENDIF 

IF ( ( QZ  .GT.  .3261}  .AND.  (QZ  .LB.  84.175))THEN 

VEL*  (  VISC/ (  DENS*2.  *RP )  )*EXP( -3  -  18657-*-.  992696*  YY- 

A 

1  .531  93E-3*YY**2~9.87059E-4*YY**3-5.78878E-4»TY**4 

♦8.551759E-5*YY**5-3.27815E-6*YY**6) 

GO  TO  1100 
ENDIF 

IF( ( QZ  .OT.  84.175)  .AND.  ( QZ . LT . 1 40 ) ) THEN 

VEL»VISC*QZ*(4  .  166667E-2-2. 3  36  3E-4*QZ4-2 .01  54E-6*Q 

Z**2 

A  -6.9105E-9*QZ**3)/(DENS*2.*RP) 

GO  TO  1100 
ENDIF 

IF ( ( QZ  .GE.  140)  .AND.  (QZ  .LT.  4.5E7))THEN 

LOGREY  *-  1  .  295  364-.  986*  YY/ 2. 303- .  0466  77*  (  YY/2 . 30  3  )  * 

•2 

4..OOI  1235*(YY/2.303)**3 
REYN0L«10.**L0GREY 
VEL«REYN0L*VISC/(2.*DENS*RP) 

ENDIF 

C  KNUDSEN  SLIP  CORRECTION 

1100  IF(RP.GT. 100.*MFP)  THEN 

EXPFACaO  . 

ELSE 

EXPFAC*EXP(-. 656*2. *RP/ MFP) 

ENDIF 

VELsVEL*  (  1  ,4-(  1 .64  4  4..552*EXPFAC)»MFP/ 

A  (  2  .  *RP ) ) 

RETURN 
END 


SUBROUTINE  ATMOS ( ZB , DENS , VISC , MFP ) 

C  SUBROUTINE  CALCULATES  ATMOSPHERIC  TEMPERATURE,  PRESSURE  AND 

VISCOSITY 

C  POR  ALTITUDES  UP  TO  84,852  METERS  (MRS  UNITS  EMPLOYED) 

REAL  LR,NUMDBN,MFP 
IP ( ZB  .LB.  11000. )THBN 
LR— .006545 
TR-288 . 1 5 
PR-101300. 

TB«TK+LR»ZB 

PB«PR»(TR/TB)»»( .034164/LK) 

00  TO  1240 
ENDIP 

IP ( ( ZB  .QT.  11000.)  .AND.  (ZB  .LE.  20000. ))THEN 
ZR» 1 1000. 

TR-216.65 

PR-22690. 

LR-0 . 

TB-TR 

PB*PR»EXP(-.034164»(ZB-ZR)/TR) 

00  TO  1240 
ENDIP 

IP ( ( ZB  .OT.  20000.)  .AND.  (ZB  .LB.  32000. ))THEN 
ZR'-20000  . 

TR-216.65 
PR-5528 . 

LR-.001 

TB-TK+LR#(ZB-ZR) 

PB-PR*(TR/TB)**( .034164/LR) 

GO  TO  1240 
ENDIP 

IP ( ( ZB  .OT.  32000.)  .AND.  (ZB  .LE.  47000. ))THEN 
ZR-32000. 

TR-228.65 
PK-888.8 
LR- . 0028 

TB-TRfLR«*(  ZB-ZR) 

PB«PR*(TK/TB)t#( . 034164/LR) 

00  TO  1240 
ENDIP 

IF ( ( ZB  .OT.  47000.)  .AND.  (ZB  .LE.  51000. ))THEN 
ZR-47000  . 

TR-270 . 65 
PR- 110. 873 
LR-0  . 

TB-TK 

P8aPR»EXP(-.034l64»(ZB-ZK)/TR) 

00  TO  1240 
ENDIP 

IF ( ( ZB  .GT.  51000.)  .AND.  (ZB  .LE.  71000. ))THEN 
ZK-51000 . 

TR-270 . 65 
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PK»66 .9218 
LK»- . 0028 
TB»TK+LKf ( ZB-ZK ) 

PB «PK*(TK/TB)##(. 034164/ LK) 

GO  TO  1240 
BHDIF 

IP( (ZB  .QT.  71000.)  .AND.  (ZB  .LE.  84852. ))THEN 
ZK«71000. 

TK»2 1 4 . 65 
PK«3. 95536 
LK»- .002 

TB*TK+LK#(ZB-ZK) 

PB»PK«(TK/TB)«»( .034164/LK) 

00  TO  1240 
BNDIP 

IP( ( ZB  .OT.  84852 ) ) THEN 
WRITS  (6,1230) 

00  TO  1241 
BNDIP 

1230  PORMAT( *  PARTICLE  ALTITUDE  BEYOND  RANGE  OF  PROGRAM 
APPLICABILITY') 

1240  VISC*1 .46E-6»TB»*1 .5/(TB+1 10.4) 

DENS*.003484»PB/TB 
NUMDEN*DENS»2.55B25/1 .23 

C  PER  CUBIC  METER:  2.55825  MOLECULES  PER  CUBE  AT  SEA 

LEVEL. 

C  DENSITY  AT  SEA  LEVEL  IS  1.23  KG/CUBE 

MFPa 1/(1 . 4 1 4*NUMDEN*4 . 5B- 1 9 ) 

C  MFP  a  1 / ( R00T2*N*SIGMA ) 

1241  RETURN 
END 

END  OF  FILE 


APPENDIX  B 


THE  LOGNORMAL  DISTRIBUTION  AND  ITS  MOMENTS 


The  lognormal  distribution  derives  from  the  normal  distribution  if  the  variable 
of  integration  is  replaced  by  its  logarithm:  if 

r  -*>  ln(r), 


Xt' 

dr  „  f-ro 

~r— — exp  —14  - 

v2n  a  <7 

dr  f/»(rH»(r«)f 

^rxprl — ? — JJ 


(B.l) 


(B.2) 


where  rm  is  the  median  radius  of  the  number  distribution  (50%  of  particles  are 
smaller  than  rm )  and  0  governs  the  dispersion  of  the  distribution: 


a  rM.l3%  rm  [  „  1 

0  »»  exp  -  —  exp  -  —  exp  S 

rm  r  15.87%  1  J 


(3.3) 


S  is  known  as  the  "logarithmic  slope"  (or  sometimes  just  "slope")  of  the  distribution. 
The  distribution  function  appears  in  several  forms: 


dN 

d{lnr) 


V.  /n(r)-/n(rm)  f 

- 1 - 

.  ~  /  J 


(B.4) 


1  j,  Mr)-/n(rm) 

y— —  exp  —n  - ^ - 

72n0r  (I  0 


(B.5) 


1  „  /n(r)-/n(rm) 

n(r)- v^rxp  -* — 0 — 


(B.6) 
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The  last  expression  is  used  in  the  present  research.  The  symbol  "  *  "  denotes  the 
function  as  normalized,  i.e. 

00 

fn(r)dr  ■»  1  (B.7) 


The  mode  radius  of  a  lognormal  distribution  lies  below  rm : 


—0* 

rmod*  m  rm  e 


The  mean  radius  lies  above  r„ 


T  -r  C*? 
’mean  1  m  c 


(B.8) 


(B.9) 


Any  moment  of  a  lognormal  distribution  is  itself  lognormal  (true  for  both  positive 
and  negative  moments).  This  can  be  demonstrated  by  successive  substitutions  for  the 
variable  of  integration.  Consider  the  l  th  moment: 


A  (r )  m  frln(r)  dr 
0 

r 


/ 


-j  i 


exp 


— % 

ln(r)-ln{rm ) 

71 

P  \ 

wf 


(B.10) 


letting 


fn(r)-/»(rm) 


P 


(B.ll) 


A(r)  becomes: 


lnr—lnrm 

1 


A{r)  -  J 


—00 


(B.12) 


completing  the  square  in  the  exponent  of  the  integrand  yields: 

lnr—lnrm 
1  “ 


A(r)-^—e‘  W2e'<'”r->  / 


(B.13) 


changing  variables  again,  let 


W  »  Z  —  1 0 


(B.  14) 
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Then: 


Inr-^nu+W*) 

X{r)  —  [e*iv  /  e-»“Vw 

Letting 

»■/  -  In  rm  +  10* 

yields  a  very  useful  integral  expression  for  A(r): 


(B.15) 


(B.16) 


M 


0-[ 


e»ivv(,nr«) 


lnr—lnri 
0 


J—  /  e-“*,dw 


(B.17) 


The  expression  above  is  readily  evaluated  using  approximations  in  Abramowitz  and 
Stegun  (72:  932-933).  Equation  B.17  can  be  manipulated  further  into  lognormal 
form  by  letting: 

lnr  —  In  rt 


w 


(B.18) 


Then, 


-Vi 


A(r)  -  e 


WV  +  l{lnrm)jC_ 


In  r'  -  lnr/ 


V27T/?/ 


-dr' 


(B.19) 


Differentiating, 


dA(r)  _  vi/a^  +  /(/„  rm 


-Vi 


la  r  -  In  r, 


-dr 


(B.20) 


dr  ~  V27T 0r 

o(r)  ( analogous  to  cqn.  B.  6) 

Thus  the  / th  moment  is  lognormal  with  the  same  dispersion  (/?)  but  a  different 
median  radius,  r( : 

r,  (B.21) 

Normalized  expressions  for  a(r)  are: 

lnr  -  lnr,  [ 

.  *  ] 


-Vi 


a(r )  »  — — — e 

V2# 


(B.22a) 
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-Vi 


'2i r/?r 


(B.22b) 


Note  that  the  average  value  of  ri  is: 


00 


<  rl  >  —  Jrln(r)dr  »  A  (oo)  «  t 


Vi /***+/  lnrm 


(B.23) 


The  ease  of  manipulating  the  moments  of  the  lognormal  distribution  greatly 
simplifies  n(r)  parameter  analysis. 


APPENDIX  C 


NATHANS'  DISTRIBUTION  AND  ITS  MOMENTS 


Nathans  fit  his  particle  size  data  to  a  function  of  the  following  form: 

exP  ,  r  <  r0 

-nB(r)-nc  ,  r>r0  (< 


(C.l) 


where  n0  «■  n(r0). 

The  transition  rsdius,  r0 ,  is  found  by  requiring  that  the  number  size  distribution 
and  its  first  derivative  be  continuous  there.  This  radius  is  given  by: 


(C.2) 


Note  that  the  3rd  moment  (mass  distribution)  of  n(r)  is  not  analytic  for  p  <  4,  i.e., 


»  fr  y 

■  To 

n‘  T 

» 


r3dr— *oo  ,  p  <4 


(C.3) 


This  problem  is  mitigated  by  defining  a  maximum  radius,  rmax  as  a  function  of  the 
mass  fraction  (/m)  below  r0.  Starting  with: 


fnA(r)r3dr 

0 _ 

r§  Turn* 

JnA  (r  )r3dr  +  /  nB(r)r3dr 


(CA) 


Solving  for  rmax: 


[r0(l  -l/fm)(p  -4)e*nr-4rV2Z0CNF[0(n-4)}  +  1  ,  p  <4 

ro  exP  [Vj&rfll  “  I m  film  ]  »  p  “  4  (c’5) 
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where  CNF  designates  the  normal  or  Gaussian  probability  function: 

CNF (*)  -  /  -i 


(C.6) 


Having  defined  rmaT  for  the  3rd  moment,  the  lower  moments  must  be  normalized 
assuming  no  particles  exist  with  radius  greater  than  rmax. 

The  normalized  moments  of  n(r)  are: 


Jnt(r)rldr 
0 

1a  +4 


Mr) 


,  r  <  r0 


JnA{r)fldr  +  fnB(r)fldr 


r. 


1a  +  h 


.  r  >r0 


where 


I  a  -fnArldr  -V5 faWr.WrPftCNF 


In M  ~  In (rt) 


/? 


and 


h  -7/  (tJ*  1  +>-1 


Tg  In 


snax 


Thus, 


y/S^e^r^r^CNF 


,  i  1  “  P  _1 


In(r)-lmi) 


A(r) 


h  +  h 


-  .  r  <  t. 


t  ■  r°  L<-r+i  .rt-p+i1 
*  .  L~P-±L-L— - : - I,  r>r0,l+p- 1 


1a  + 


(C.7a) 


(C.7b) 


(C.8) 


(C.fla) 

(C.Ob) 


(C.lOa) 


(C.lOb) 


I  A  +  r$  In  (r/r0) 
(4  +  '* 


,  r  >  t0  ,  /  —  p  —  1 
140 


(C.lOc) 


tatifitflium 


r*  ’\  -  v  ■*  ,  *  *  < 

-  -These  equations  are  the  basis  for  the  activity  group  computation  in  subroutine 
'SIZNAT  (see  Appendix  E  for  a  listing). 
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APPENDIX  D 
Program  SHOT60  Listing 
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PROGRAM 

SH0T60( INPUT, OUTPUT, TAPE5 a IN PUT ,TAPE6*OUTPUT, TAPE2 ,TAPE4 , 
*PLFILB*0 ) 

C  PROGRAM  IDENTIFIES  BEST  VALUES  OF  MOMENT,  MEDIAN  RADIUS, 

AND  BETA  OF 

C  "LOG  NORMAL  STRATOSPHERIC  PARTICLE  SIZE  DISTRIBUTION  BASED 
UPON 

C  "ATMOSPHERIC  TEST  DATA 

DIMENSION  TRUED A(20) ,BURDEN(20) ,RC0UNT(336) , MONTH (336) 
DIMENSION  R0(20) ,R(2C) ,TRUEYR(20) ,YEAR(336) 

DIMENSION  A(8  v8) ,B( 8 , 8) ,DELTA(8,8) 

COMMON  RGROUP ( 60 ) , WEIGHT (60) .OFFSET ( 1 1 ) , Y ( 1 1 ) ,MC( 11) 

REAL  MO,  M0 1 ,  MONTH,  MC ,  LOGR 
INTEGER  OPPSET 
C  OPEN  AUXILIARY  DATA  FILE 
OPEN ( 2 , FILE* ' PI ’ ) 

C  FIRST  ESTABLISH  VALUES  POR  DATA  (TRUE)  AGAINST  WHICH  TO 

OPTIMIZE 

TRUED A ( 1 ) *1 .60 
TRUED A ( 2 )  *  1 .20 
TRUEDAC 3) * -903 
TRUED A ( 4 )  *  1 .00 
TRUEDAC5) *.933 
TRUED A ( 6 ) *  1 .62 
TRUED A ( 7 ) *  1  .10 
TRUEDAC  8) * . 933 
TRUED A ( 9 ) *2 . 07 
TRUEDAC  1 6 ) *6 . 4  1 
TRUEDAC 1 1 ) *3 .56 
TRUEDAC 12) *1 . 47 
TRUEDAC 13) » -750 
TRUEDAC 1 4 ) » . 3 1 0 
TRUEDAC 15)*. 280 
TRUEYRC 1 )«6./12. 

TRUEYR(2)*12./12. 

TRUEYR(3)*24./12. 

TRUEYRC4)*36./12. 

TRUEYRC5)*48./12. 

TRUEYR(6)*60./12. 

TRUEYR(7)*72./12. 

TRUEYRC8) *84. /12 . 

TRUEYR(9)*96./12. 

TRUEYRC 10)*108./12. 

TRUEYRC 11)*120./12. 

TRUEYRC 12 ) *  1 32 . / 12  . 

TRUEYRC 13)*144./12. 

TRUEYRC 14)  =  156  .  /  1  2  . 

TRUEYRC 1 5 ) *  1 68 . /  12  . 

C  INPUT  OFFSET, YIELD, SR90  FOR  ATMOSPHERIC  "EVENTS” 

OFFSETC 1 ) *6 
0FFSETC2) *58 
OFFSETC  3) *64 
OFFSETC  4) *82 
0FFSETC5) *106 


OPFSBT  <  6 ) *  1 14 
OPFSBT( 7 ) ■  1  86 
OPFSBT( 8 )  •  1  88 
OFFSET ( 9 ) *204 
OFFSBT( 10)*210 
OPFSST( 1 1 ) *2 1 0 
Y<1)«9.6 
tY(2)«3.5 
T ( 3 ) *  1 .5 
Y( 4 ) *2 . 4 
Y( 5) *  .  7 
1(6) *3*6 
Y ( 7 ) *2 . 8 
Y( 8 ) *4 1 .5 
Y ( 9 ) *5 ♦ 3 
Y(10)*25. 

Y( 1 1 ) «6 . 8 
MC( 0*2.4 
MC( 2 ) « . 88 

MC ( 3 ) * • 3 
MC ( 4 ) *  1 . 4 
MC ( 5 ) *  1  .  1 
MC( 6) *  1  .  1 
MC ( 7 ) * . 644 
MC( 8 ) * . 689 
MC(9)a.427 
MC ( 1 0 ) *4 . 56 
MC ( 1 1 )«1 .95 

SET  VALUES  FOR  THE  MOMENT,  MEDIAN  RADIUS,  AND  SLOPE  (BET) 
BETA  IS  THE  NATURAL  LOO  OF  THE  SLOPE 

WRITE ( 6  ,  • ) '  INPUT  NDIST  (1*L0G  NORM,  2=NATHANS  HYBRID)' 
READ ( 5 , • ) NDIST 
WRITE ( 6 , • ) '  INPUT  MOMENT' 

READ ( 5 , • ) M0 1 
IsO 

DO  20  BEIal .01 ,4.  ,  .2 
1*1  +  1 
J*0 

DO  10  LOORa-2 .  ,  1  .  ,  .  1 
J«J+1 

RMIalO .••LOOR 
FORMAT  (A) 

IP(NDIST.EQ.2)WRITE(6  ,•) '  INPUT  FMASS  (P1*4)  ' 
IF(NDIST.EQ.2)READ(5  ,  • )  F 
F0«0. 

CALL  STRATF (NDIST, MO  1 , RM 1  ,BE1 , F , BURDEN , RCOUNT ) 

DO  30  Ks  1  0  ,  1  5 
RO(K)*BURDEN(K)-TRUEDA(K) 

RO(K)*RQ(K)#RO(K) 

FO*FO+RO(K) 

CONTINUE 

WRITE ( 6 , 32 )  I , J  ,  BE  1  ,RM1  ,F0 
WRITE(2,50)I,J,L0GR,BE1 ,F0 
FORMAT (  I4,I4,E12.5,E12.5,E12.5) 


10  CONTINUE 

20  CONTINUE  „  „ 

32  FORMAT  (  •  I '  ,  14  ,  '  J',I4,'  B',E8.3,’  R’,E8.3,’  FQ’,E8.3) 

DO  40  I«1  ,336 
MQNTH(I)*PL0AT(I)/2. 

TEAR ( I ) *MONTH{ I )  /  1  2  . 

40  CONTINUE 
00  TO  999 

C  DISSPLA  PLOTTING  CALLS  FOLLOW 

77  CALL  COMPRS 

CALL  XNAME('TEAR$' ,100) 

CALL  TNAME( 'STRATOSPHERIC  BURDEN  (MC)$’,100) 

CALL  PAOE ( 8.5, 11.) 

CALL  AREA2D ( 5.1,4.) 

CALL  HEADINCSR90  BURDEN  FROM  U.S.  FOREIGN 
TESTS*’ , 100, 1 . , 1 ) 

CALL  INTAXS 

CALL  QRAF(0 . , 1 . , 14. , 0, , 1 . ,8. ) 

CALL  THKFRM( .02) 

CALL  FRAME 
CALL  MARKER ( 15) 

CALL  CURVE(TRUEYR,TRUEDA,15, 1 ) 

CALL  MARKER ( 0 ) 

CALL  CURVE( TEAR.RCOUNT, 336, 12) 

CALL  GRID ( 1,1) 

ENCODE (60, 1 5, LABEL) M0 1 , RM1 ,BE1 ,F0 
15  PORMAT( 'MOMENT*'  ,F4 . 1 , *  RM*',F4.1,»  SLOPE *',  F4 . 1  , 

*  '  FOa'  ,89.3,'$') 

CALL  RLMESS(LABEL, 100, 1 . ,7. ) 

CALL  XDTAXS (540101 , 4HTEAR , 680 1 0 1  ,5.1, '$’,100,0. ,-.75) 

CALL  ENDPL(O) 

CALL  DONEPL 
999  END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCC 

SUBROUTINE  MTXEQ ( A , X , B , N , K ) 

C 

C  MATRIX  EQUATION  SOLVER 

C 

C  USAGE... 

C 

C  TO  SOLVE  THE  LINEAR  SYSTEM  AX  =  B 

C 

C  CALL  MTXEQ ( A , X , B , N , K ) 

C 

C  WHERE  A  MUST  BE  DIMENSIONED  N  X  N 

C  X  MUST  BE  DIMENSIONED  N  X  K 

C  B  MUST  BE  DIMENSIONED  N  X  K 

C  N  IS  THE  NO.  OF  EQUATIONS  (ROWS  IN  A,X,B) 
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ooo  o  o  o  o  o  o  o  o  o  a  a  a  a  a  a  oooonoorioo 


K  IS  THE  NO.  OF  SOLUTION  VECTORS  (COLS.  IN  X,B) 


NOTE...  TO  CHANGE  DIMENSIONS  OP  ARRAYS  C  AND  PIV,  ALSO 
CHANGE  VALUES  OF  NMAX  AND  NRMAX  IN  DATA  STATEMENT. 


DECLARATIONS 

INTEGER  RPjNP,J,I,IFR0M,IP1,IPIV,IT0,L,NP1 ,NPJ,NPK 
REAL  A(8, 8) ,B(8, 8) ,X(8,8) 

REAL  PIV( 16) ,C(8, 16) , ATPB, RM 
DATA  NMAX,  NKMAX/  8,  16/ 

GET  ARGUMENTS  N  AND  X 

NP»N 

KP»K 

TEST  N  AND  K  FOR  CORRECT  RANGE 

IF(NP.LE.O.OR.NP.GT.NMAX)  GO  TO  190 
IFUP.LE.O.OR.  (NP+KP)  .GT.NKMAX)  GO  TO  190 

MOVE  ARRAYS  A(I,J)  AND  B(I,J)  INTO  C(I,J) 

DO  10  Ja 1 ,NP 
DO  10  I  a  1  , NP 
10  C(I, J)aA(I, J) 

DO  20  Jal ,KP 
NPJaNP+J 
DO  20  Ial ,NP 
20  C(I,NPJ)aB(I, J) 

SET  TO  PERFORM  N  ELIMINATION  SWEEPS  (Ia1,N) 

NP1aNP+1 
NPKaNP+KP 
DO  120  Ial ,NP 
IP1aI+1 

SEARCH  FOR  NEXT  PIVOT  ROW  (I-TH  PIVOT  IS  IN  COL.  I) 

ATPEaO. 

DO  40  Jal.NP 

IF  (ABS(C(J,I) )«ATPE)  40,30,30 
30  ATPEa ABS ( C( J  ,  I ) ) 

IPIVsJ 

40  CONTINUE 

OPERATE  ON  THE  PIVOT  ROW 


IF  ( ATPE )  210,210,50 


non  non  non  o  o  o 


50 


DO  60  J«IP1 , NPK 

IF  (CCIPIV.I)  .EQ.  0)  THEN 
WRITB( 6 | 302 ) 

WRITE(6 ,327) 

STOP 
ELSE 

PIV(J)aC(IPIV, J)/C(IPIV,I) 

BNDIF 

60  CONTINUE 

PERFORM  ELIMINATIONS  BELOW  THE  DIAGONAL  (COL.  I) 

IFROMaNP 
ITO*  NP 

0  IP  (IFROM-IPIV)  80,100,80 
0  RM»-C( IFROM  ,  I ) 

DO  90  J.IP1 , NPK 
90  C(ITO,J)aC(IFROM,J ) ♦RM*PI V ( J  ) 

ITOalTO- 1 

100  IFROMalFROM- 1 

IF  ( IFROM-I )  110,70,70 

POT  THE  I-TH  PIVOT  ROW  IN  THE  VACATED  ROW  I 

110  DO  120  J.IP 1 , NPK 
120  C( I , J) aPIVC  J) 

NOW  DO  THE  BACK  SOLUTION 

I  a  NP 

130  IP1*I 

IaI-1 

IF  (I)  160,160,140 
140  DO  150  JaNPI , NPK 
DO  150  L«IP1 , NP 

150  C(I,J)«C(I, J)-C(I,L)»C(L,J> 

GO  TO  130 

MOVE  THE  SOLUTION  TO  ARRAY  X(I,J) 

160  DO  170  Ja 1 , KP 
NPJaNP+J 
DO  170  I  a  1 , NP 
170  X ( I , J ) aC ( I , NP J ) 

180  RETURN 
C 

190  WRITE{6, 1000)  NP,  KP 
C  CALL  SYSTEM  (  200,  'L’  ) 

STOP 

210  WRITE(6, 1001  ) 

C  CALL  SYSTEM  (  200,  'L’  ) 

STOP 

302  FORMAT ( 4 ( / ) , 5 1 H  DIVISION  BY  ZERO  OCCURED  IN  SUBROUTINE 

MTXEQ  - 
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•,32H  DURING  EVALUATION  OF  PIV(J)  •••) 

327  FORMAT( 24h  PLEASE  CHBCK  YOUR  INPUT, 4(/)) 

1000  FORMAT (3HON*,I12,3H  K*,I12,35H  ARE  INCORRECT  FOR  SUBROUTINE 
MTBQ  ) 

1001  FORMAT  ( 37H0DBT ( A ) «0  IN  CALL  TO  SUBROUTINE  MTXEQ) 

END 

C  SUBROUTINE  SYSTEM 

C  RETURN 

C  END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCC 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

CCCCCCCCCCCCC 

SUBROUTINE 

STRATF (NDIST, MOMENT, RMICRON, SLOPE, F .BURDEN, RCOUNT) 

C  VERSION  TO  GENERATE  PSEUDO  DATA  TO  CHECK  OPTIMIZATION  TECHN 

C  TAPES  FROM  KEYBOARD , TAPES  TO  SCREEN , TAPE2»EVENT 

DATA ,TAPE3»DATA  OUT  FOR  PBBBR 
C  PLOT 

DIMENSION  CONTR ( 60 ) ,  ZGR0UP(60) 

DIMENSION  GRPSUM( 1 1 , 336) ,  RCOUNT(336),  GCOUNT(336), 

FALL ( 336 ) 

DIMENSION  PLTEAU(336),  TVEL(60) 

DIMENSION  BURDEN(20) , HOC (  1  1) 

COMMON  RGROUP ( 60 ) .WEIGHT (60) , OFFSET ( 1 1 ) , Y ( 1 1 ) , MC ( 1 1 ) 

REAL  MONTH, LNY, MOMENT, MC 

INTEGER  OFFSET, EVENT, CONTR, TABS, TAB 

EXTERNAL  STRAP, VEL 

RM*RMICR0N»1 .E-6 

B= ALOG ( SLOPE ) 

C  OPEN  AUXILIARY  DATA  FILES 

C  OPEN (2 ,FILE»» EV60' , STATUS a » OLD »  ) 

C  0PEN(3,FILEs»PSEUDA* ) 

C  INPUT  RUN  PARAMETERS  FROM  KEYBOARD 

C  WRITE(6,5)  *  WHAT  IS  TROPOPAUSE  ALT  (M)  ?• 

ZTROP* 12000 . 

C  READ ( 5 , • )  ZTROP 

C  WRITE(6,6)ZTROP 

5  FORMAT ( A ) 

6  FORMAT (  *  TROPOPAUSE  ALTITUDE  (M)**fE10.3) 

C  INITIALIZE  ARRAYS 

DO  10  TABSsI ,336 

PLTEAU(TABS)sO. 

RCOUNT { TABS ) *0 . 

GCOUNT ( TABS ) sO . 

FALL ( TABS ) »0 . 

10  CONTINUE 

DO  20  EVENTsI ,  1 1 
DO  30  INCRE* 1 ,336 

30  GRPSUM(EVENT(INCRE)=0. 

20  CONTINUE 

C  DEFINE  PHYSICAL  CONSTANTS 

G*9 . 8 

C  SET  UP  SIZE  GROUPS,  LOG  NORMAL,  OR  LOG  NORMAL  +  EXPONENTIAL 
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TAIL 


HGR0UP»60 

IF  ( NDIST . EQ . 2 )CALL  SIZHA25 ( RM , B , MOMENT , F ) 

IF( NDIST . EQ . 1 ) CALL  SIZCNF ( RM , B , MOMENT ) 

C  INPUT  TERMINAL  VELOCITIES  FOR  EACH  GROUP  AT  TROPOPAUSE 

DO  190  J*1 ,  NGROUP 
TVEL( J )*VEL(ZTROP,RGROUP(J)) 

C  WRITE( 6 » 1 9 1 )  RGROUP( J) ,TVEL( J) 

190  CONTINUE 

191  FORMAT(  2E10.3) 

C  DEFINE  DETONATION  PARAMETERS 

NBURSTa 1 1 
HOC ( 1 ) *26000  . 

HOC ( 2 ) s2 1 000  . 

HOC ( 3 ) *20000 . 

HOC ( 4 ) a2 1 000 . 

HOC ( 5 ) *  1 7000 . 

HOC (6)»23000. 

HOC (7 ) *22000 . 

HOC ( 8 ) a  38000 . 

HOC ( 9 ) *25000  . 

HOC (  1 0 ) *35000  . 

HOC  C 1 1 ) *26000  . 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
C  FOR  EACH  EVENT  COMPUTE  STRATOSPHERIC  BURDEN  VS  TIME 

DO  710  EVENT*7 , 1 1 

C  COMPUTE  STABILIZATION  ALTITUDE  AND  INITIAL  Z 

DISPERSION.  Y  IN  MT 

C  USING  HOPKINS  FORMULA 

LNYaALOO(Y( EVENT)* 1000.  ) 


C 

Z0MaEXP(7 -889+. 34*LNY+.001226»LNY»#2-. 005227»LNY»»3 
C  *  +.000417#LNY*#4) 

Z0M=H0C( EVENT) 

SIGZ0*2 . 15E3 

C  DETERMINE  WHICH  SIZE  GROUPS  HAVE  HIT  THE  GROUND  WITHIN 

1  DELTAT 


DELTAT*365 .25/24 . 
DELSECaDELTAT»24 ,»3600 . 

DO  420  J  a  1 f NGROUP 
TFALMXaZOM/TVEL( J) 
IF(TFALMX.LT.DELSEC)THEN 
CONTR( J)aO 
ELSE 

CONTR ( J ) a  1 


ENDIF 

I F ( (CONTR(J) . EQ.O)  .AND.  ( CONTR ( J- 1 ). EQ . 1 )) THEN 


JCUTaJ-1 

ENDIF 

420  CONTINUE 

C  SET  INITIAL  ALTITUDE  OF  EACH  GROUP  TO  ZOM 

DO  440  J  a  1  , NGROUP 
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ZGROUP ( J ) aZOM 
440  CONTINUE 

C  FIND  INITIAL  FRACTION  OF  EACH  CONTRIBUTING  GROUP  IN 

THE  STRATOSPHERE 


C 

479 


HT  (  J )  • 


INCREa 1 
TLAPSBbO. 

TABSbOFFSBT(EVENT) 

VRITB(6 ,479)  ZGROUP(I) , SIGZO , TLAPSE 
FORMAT (  'ZGROUP, SIGZO, TLAPSE’ ,3E10.3> 

FRACaSTRAF ( ZGROUP ( 1 ) , SIGZO , TLAPSE , ZTROP ) 

DO  480  J b 1 ,  JCUT 

GRPSUM ( EVENT, INCRE) *GRPSUM( EVENT, INCRE )+WEIG 


480 


495 

C 

MC ( EVENT )  , 
C 


FRAC*MC (EVENT) 

CONTINUE 

RCOUNT( TABS) aRCOUNT ( TABS )+GRPSUM( EVENT, INCRE) 
DO  495  TABbOFFSET(EVENT) ,  336 

PLTEAU(TAB) aPLTEAU( TAB) +MC( EVENT )- 
GRPSUM(EVENT, INCRE) 

CONTINUE 

VRITSC  6 , 500 )  FRAC,  TABS,  RCOUNT ( TABS ) , 

GRPSUM ( EVENT , INCRE ) ,PLTEAU( TABS) 


500  FORMAT  ( 

•FRAC' ,B 10. 3, 'TAB' , 14 , 'COUNT' , E 1 0 . 3 , ' MC  ’  , 

*  E10.3,'GR0UPSUM'  , El  0.3, 'PLATEAU'  , El  0-3) 

C  SET  TIME,  INCREMENT  IN  STEPS  OF  1/24  YEAR  (TWO  WEEKS) 

510  INCREsINCRE+l 

TABS  s  IN  CRB<f  OFFSET  (EVENT)-I 
IF(TABS.GT.336)THEN 
GO  TO  710 
ENDIF 

IFC (GRPSUM (EVENT, INCRE- 1 ) ) .LT. 1 . E- 1 ) THEN 
GO  TO  710 
ENDIF 

TLAPSEaINCRE*DELTAT 

C  COMPUTE  NEW  Z(T)  FOR  EACH  CONTRIBUTING  GROUP 

DO  680  J  a  1 , JCUT 

RPsROROUP(J) 

IF(ZQROUP(J) .GT.O. ) THEN 
ZBsZGROUP ( J ) 

ELSE  ZBbO. 

ENDIF 

DZbVEL(ZB,RP)»DELSEC 

ZGROUP(J)sZGROUP(J)-DZ 

GR PSUM ( E VENT, INCRE) sGR PS UM( EVENT , INCRE) f 


STRAF ( ZGROUP ( J ) , SIGZO , TLAPSE , ZTROP ) •WEIGHT ( J ) 

A  *MC ( EVENT ) 

C  WRITE(6,679)EVENT, INCRE, J, ZB, VEL(ZB,RP) , 

C  *  STRAF( ZGROUP( J) , SIGZO, TLAPSE, ZTROP) , 

C  *  GRPSUM( EVENT , INCRE) 

679  FORMAT (  '  EV',13,'  INCRE’, 14,’  J',I3, 

*  ’  ZB ’ , E 1 0 . 3  »  '  V EL ' , E 1 0 . 3  ,  '  STR AF ’ , E 1 0 . 3  ,  ' 
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GSUM ’ , BIO . 3) 

680  CONTINUE 

RCOUNT( TABS )*RCOUNT( TABS )+GRPSUM( EVENT , INCRE)* 
*  SXP(-6.78B-5»TLAPSE) 

700  00  TO  510 

710  CONTINUE 

DO  810  1*1 , 15 
K*I»24-24 

BURDBN(I)*RCOUNT(K) 

810  CONTINUE 

850  FORMAT (  P5.1 ,E10. 3,810.3) 

C  ENDFILE  3 

RETURN 
END 


SUBROUTINE  SI 2CNF ( RM , 8 , MOMENT ) 

C  PROGRAM  TO  PARTITION  PARTICLE  MASS  DISTRIBUTIONS  INTO 

GROUPS 

COMMON  RCBNTR ( 60 ) ,MASPRA(cC ) ,OFFSBT( 11), 1(11), MC (11) 
DIMENSION  RLEFT ( 60 ) ,RRIGHT(60) ,MASCNF(60) 

REAL  MASCNF, MASFRA, MOMENT, MC 
INTEGER  GROUP, OFFSET 
EXTERNAL  CNF 

C  CNF  IS  THE  CUMULATIVE  NORMAL 

FUNCTION (RADIUS.RM, BET A, MOMENT) 

C  0PEN(7 ,FILE«* TAPE7’ ) 

C  DEFINE  LOG  NORMAL  FUNCTION  PARAMETERS,  RM  BETA 

30  FORMAT  (13 ,4(E10 . 3) ) 

C  WRITE  (6,40) 

40  FORMAT ( '  GROUP  RLEFT  RRIGHT  RCENTROID  MASS  FRAC') 

EXPON  =  -9 . 

GROUPal 

RLEFT (GROUP ) =0 . 

RRIGHT(GR0UP)s10. ••EXPON 
RCENTR(GROUP) a RRIGHT (GROUP) 

C  WRITE ( 6,60)  RRIGHT (GROUP) ,RM,B 
C60  FORMAT  (  3E10.3) 

MASCNF( GROUP) aCNF( RRIGHT (GROUP) , RM , B , MOMENT ) 

MASFRA ( GROUP ) *MA3CNF( GROUP) 

DO  130  GR0UP=2,60 

£XPONaEXPQN+. 10 

RRIGHT ( GROUP) =10 .••EXPON 

RLEFT ( G ROUP ) = RRIGHT (GROUP- 1 ) 

RCENTR( GROUP) *SQRT( RLEFT (GROUP) ‘RRIGHT (GROUP ) ) 
RsRRIGHT(GROUP) 

MASCNF( GROUP ) aCNF(R , RM,B , MOMENT) 

MASFRA ( GROUP ) aMASCNF( GROUP ) -MASCNF (GROUP- 1 ) 

C  WRITE  (6,30) 


151 


ORO0P,RLEFT(GROUP) ,  RR IGHT ( GROUP ) , RCENTR ( GROUP ) 
C  *  ,MASFRA ( GROUP ) 

130  CONTINUE 
RETURN 
END 

ClltllltltMMIMMtlllMMtMIMIIIItllltMft 
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•!*•*•••••••••• 

SUBROUTINE  SIZNA25 ( RH , B , MOMENT , F ) 

COMMON  RCBNTR(60) ,ACTFRA(60) 

C  SUBROUTINE  PARTITIONS  NATHANS  R4  DISTRIBUTION  INTO  ACTIVITY 

GROUPS 

DIMENSION  RLEFT ( 60 ) .RRIGHT (60) ,ACTCUM(60) 

REAL  MOMENT, PI,MC 
INTEGER  GROUP, OFFSET 
EXTERNAL  CNF 

C  CNF  IS  THE  CUMULATIVE  NORMAL  FUNCTION 

C  DEFINE  NATHANS  PARAMETERS 

R0«RM»EXP(3.*B«B) 

R25jRM»EXP(S  5*B*B ) 

PI*3. 14159 
R2PI*SQRT(2.»PI) 

RMAXaR0»EXP(R2PI»B»( 1 . -F)/2./F) 

30  FORMAT  (I3,5(E10. 3)  ) 

WRITE ( 6 , 40 ) 

40  '  FORMAT ( •  GROUP  RLEFT  RRIGHT  RCENTR  ACTV  FRAC 
CUM'  ) 

EXPON  a-9  .  1 
GROUPsO 

60  FORMAT  (  3E10.3) 

DO  280  GROUPal ,  49 
EXPONa  EXPON* .  1 
RRIGHT( GROUP) a  10 .••EXPON 
IF  (GROUP. EQ.1)  THEN 
RLEFT(GR0UP)*0. 

RCENTR ( GROUP ) a RRIGHT ( GROUP  ) 

ELSE 

RLEFT( GROUP ) a RRIGHT (GROUP- 1 ) 

RCENTR(GROUP)aSQRT(RLEFT(GROUF)»RRIGHT (GROUP)) 

END  IF 

RT3RRIGHT(GROUP ) 

C  CALCULATE  DENOMINATOR  WHICH  IS  PROPORTIONAL  TO  TOTAL 

ACTIVITY 

DENOM*R2PI»B»RM»RO»»2.5*EXP(3. 1 2 5*B»B ) • CNF ( RO , RM , B , MOM 

ENT) 

*  +R0*#4*2.*( 1 ./SQRT(R0)-1./SQRT(RMAX)) 

C  NOW  CALCULATE  NUMERATOR  WHICH  IS  PROPORTIONAL  TO  ACTIVITY 

DOWN 

IF  (RT.LT.RO)  THEN 

CNF P A RT = CNF ( RT, RM,B, MOMENT) 

ELSE 

CNFPARTsCNF(RO  ,  RM , B , MOMENT) 

END  IF 
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STaR2PI*B*RM*R0**2.5*EXP( 3- 125*B*B)*CNFPART 
BTaR0**4*2.*( 1 ./SQRT(R0)-1 ./SQRT(RT) ) 

IF  (RT.LT.RO)  BT*0. 

XNUMsST+BT 

ACTCUH(GROUP) aXNUM/DENOM 
IF  (OROUP. EQ. 1 )  THEN 

AC TFRA(GROUP)aACTCUM (GROUP) 

ELSE 

ACTFRA(GROUP)aACTCUM(OROUP)-ACTCUM (GROUP-1) 
END  IP 
VRITE( 6 , 30 ) 

GROUP, RLEFT ( OROUP ) ,RRIGHT( GROUP) , RCENTR ( GROUP ) 

A  , ACTFRA(GROUP) ,ACTCUM (GROUP) 

280  CONTINUE 

C  DEFINE  GROUP  WHICH  LUMPS  BIG  PARTICLES  TOGETHER  WITH 

RCENTRa 1 00  MICRONS 
0R0UP*50 

RLEFT (GROUP) aRRIGHT ( GROUP- 1 ) 

RRIGHT(GROUP)*RMAX 

RCENTR (GROUP) a RRIGHT( GROUP- 1 ) 

ACTPRA ( GROUP ) *  -ACTCUM( GROUP- 1  ) 

ACTCUM ( OROUP  TCUM ( GROUP- 1 ) +ACTFR A ( GROUP ) 

WRITE ( 6 , 30) 

GROUP,  RLEFT  (  GROUP)  ,  h I  JtlT(  GROUP)  ,  RCENTR  (  GROUP  ) 

*  , ACTFRA( GROUP ) , ACTCCM( GROUP ) 

RETURN 

END 

CMIIIIMIIMMMIIIIMlfHHMMIIIMMMIMMmilKtlMI 
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FUNCTION  CNF(R,RM,B, MOMENT) 

C  COMPUTES  CUMULATIVE  NORMAL  FUNCTION  FOR  MASS  DISTRIBUTION 

REAL  MOMENT 
R3sRM*EXP(M0MENT*B*B) 

Z=(ALOG(R)-ALOG(R3) )/B 
C  WRITE ( 6 , 50  )  Z 

C  50  FORMAT  ( '  Za* ,B10.3) 

IF  (Z)  1060,1050,1050 

1050  CNF=1 . - . 5 / ( 1 .♦. 196854»Z+. 1 1 5 1 9 4 • Z • * 2  ♦ . 000 344»Z»» 3 

♦ . 01 9527*Z 

A  «*4  )*«4 


GO  TO  1070 
1060  Za-Z 

CNP*  .  5/  (  1  196854*  Z+.  1 15194*  Z*«2*.  000  344  »Z»»3+.  01  9527 

*Z**4 ) 


•*  4 


1070  RETURN 
END 


C 


C 


FUNCTION  STRAF(ZB,SIGZO,TLAPSE,ZTROP) 
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REAL  KDZ 
KDZ*  .  5 

Z3CALE.6.5E3 

SIGZ»SI0Z0+SQRT(2.*KDZ»TLAPSE»24 .*3600. ) 

DELTZ«(ZB-SIGZ*SIGZ/ZSCALE)-ZTRQP 

AaDELTZ/SIGZ 

IF  (A)  1060,1050,1050 

1050  CUMNF«1 .-.5/( 1 .♦. 196854*A+. 1 15194*A**2  ♦  . 000344*A**3 
♦ . 0 1 9527*A 

*  •  •  4  )»*4 

00  TO  1070 
1060  A*-A 

C0MHF«.5/(  1  196854*  A+.  1 1  5  1  94*  A**2-(- .  000344 *A**  3-*- .  01  952 

7*A**4) 

•  •4 

1070  STRAPaCUMNF 
RETURN 
END 

CtllllMMIMIIKIMtIMIIIII.MMIIMIIiailllMIIIMMIIIIMMM 


c 


FUNCTION  VEL ( ZB , RP ) 

C  COMPUTES  TERMINAL  VELOCITY  OF  SPHERICAL  PARTICLES  USING 

FULL  DELFIC 

C  FALL  MECHANICS  AND  U.S.  STANDARD  ATMOSPHERE  EQUATIONS 

C  SUBROUTINE  ATMOS  RETURNS  DENSITY,  VISCOSITY  AND  MEAN  FREE 

PATH  AT  ALT  ZB 

REAL  MFP , LOQREY 
C  EXTERNAL  ATMOS 

CALL  ATMOS (ZB,DENS,VISC,MFP) 

DENSF»2600. 

G  a  9  •  8 

QZ«4*DENS*(D£NSF-DENS)*G*(2.»RP)**3/(3.*VISC*VISC) 

WRITE ( 6 , 1051 ) QZ, DEN3 , VISC, MFP 
1051  FORMAT ( * QZ , DENS , VISC , MFP  *,4(B10.3)) 

YYaALOO(QZ) 

C  DENSF  IS  PARTICLE  DENSITY,  G  IS  GRAV,  QZ  IS  DAVIES  NUMBER 

DEFINED  IN 

C  DNA  TR-5 1 59F-  1  ,  P24 

IF ( QZ  .LE.  .3261 )THEN 

VEL«VISC*QZ/(24.*DENS*2.*RP) 

00  TO  1100 
ENDIF 

IF ( ( QZ  .OT.  .3261)  .AND.  ( QZ  .LB.  84.175))THEN 

VEL*(VISC/(DENS*2.*RP  )  )*EXP(-3  .  1  8657-f .  99 269 6* YY- 

A 

1. 531932-3* YY**2 -9. 870592-4* YY»* 3-5. 78878E-4*YY**4 

•*-8.551759E-5*YY»»5-3. 278 1 5E-6*YY**6) 

00  TO  11 00 
ENDIF 

IF ( ( QZ  .OT.  84.175)  .AND.  ( QZ . LT . 1 40 ) ) THEN 

VEL»VISC*QZ*(4. 166667E-2-2 . 3363E-4*QZ+ 2 . 0 1 54E-6*Q 
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1**2 

*  -6.9105B-9*QZ»*3)/(DENS*2.*RP) 

00  TO  1100 
ENDIF 

IF ( ( QZ  .GE.  1  MO )  .AND.  (QZ  .LT.  4.5E7))THEN 

L0GREY«-1 .29536+. 986* YY/2.303-. 046677* (YY/2. 303)* 

•2 

‘  +.00l1235*(YY/2.303)**3 

REYN0La10.##L0GREY 
VEL«REYN0L*VISC/(2.*DENS*RP) 

ENDIF 

C  KNUDSBN  SLIP  CORRECTION 

1100  IF(RP.OT.  100.*MFP)  THEN 

EXPFACaO. 

ELSE 

EXPFAC*EXP(-. 656*2. *RP/MPP) 

ENDIF 

FEL« VEL* ( 1 .+(1 .644+.552*EXPFAC)*MFP/ 

( 2 . *RP ) ) 

RETURN 

END 


SUBROUTINE  ATMOS ( ZB , DENS , VISC , MFP ) 

SUBROUTINE  CALCULATES  ATMOSPHERIC  TEMPERATURE,  PRESSURE  AND 
ISCOSITY 

FOR  ALTITUDES  UP  TO  84,652  METERS  (MRS  UNITS  EMPLOYED) 

REAL  LK , NUMDEN , MFP 
IF ( ZB  .LE.  11000. )THEN 
LK*-. 006545 
TK»288 . 15 
PK«101 300 . 

TBaTK  +  LK*  ZB 

PB»PK*(TK/TB)**( . 034 1 64/LK ) 

00  TO  1240 
ENDIF 

IF ( ( ZB  .OT.  11000.)  .AND.  (ZB  .LE.  20000 . ) ) THEN 
ZKal 1000. 

TK«2 16.65 
PK*22690. 

LK»0 . 

TBaTK 

PB»PK»EXP(-.034164»(ZB-ZK)/TK) 

GO  TO  1240 
ENDIF 

IF ( ( ZB  .GT.  20000.)  .AND.  (ZB  .LE.  32000. ))THEN 
ZK«20000 . 

TK»2 1 6 . 65 
PK«5528  . 

LK  a . 001 

TB*TR.+LK*(ZB-ZK) 
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PB.PK»(TK/TB)«»(  .034164/LK) 

00  TO  1240 
ENDIF 

IF((ZB  .OT.  32000.)  .AND.  (ZB  .LE.  47000.  ))THEN 

.-\ZK«32000 . 

TK*228 . 65 
PK-888.8 
LEs . 0028 

TB«TK>LK* ( ZB-ZK) 

PB*PE*(TK/TB)##( . 034164/ LK) 

,0.  00  TO  1240 

ENDIF  „ 

IP  ( ( ZB  .OT.  47(500.)  .AND.  (ZB  .LE.  51000.  ))THEN 

ZK*47000. 

TK*270 . 65 
PE* 1 1 0 . 87  3 
LEsO . 

TB*TK 

PB*PK«EXP (-.0341 64«( ZB-ZK )/TK) 

QO  TO  1240 

EN01F  . 

IF((ZB  .OT.  51000.)  .AND.  (ZB  .LE.  71000.))THEN 

ZK«51000. 

TK-270.65 
PK.66.9218 
LK»- . 0028 
TB*TK+LK«(ZB-ZK) 

PB*PK»(TK/TB)»»( .034164/LK) 

00  TO  1240 
BMDXF 

IF ( (ZB  .OT.  71000.)  .AND.  (ZB  .LE.  84852. ))THEN 
ZK»7 1 000 . 

TK*214 .65 
PK. 3- 95536 
LKa- . 002 

TB*TK*LK«(ZB-ZK) 

PB*PK»(TK/TB)»»( .034164/LK) 

00  TO  1240 
ENDIF 

IF ( ( ZB  .QT.  84852)) THEN 
WHITE  (6,1230) 

00  TO  1241 
ENDIF 

1230  FORMAT ( 1  PARTICLE  ALTITUDE  BEYOND  RANGE  OF  PROGRAM 
APPLICABILITY*) 

1240  VISC.1 . 46E-6*TB## 1 . 5/ ( TB+11 0 . 4 ) 

DENS«.003484»PB/TB 

NUMDEN«DENS*2 .55E25/1  .  23  „„„  c_. 

PER  CUBIC  METER:  2.55E25  MOLECULES  PER  CUBE  AT  SEA 

E  V  EL 

DENSITY  AT  SEA  LEVEL  IS  1.23  KG/CUBE 
MFP  *  1/(1 ,414*NUMDEN*4.5E-19) 

C  KFP»1/(R00T2«N*SIGMA) 

1241  RETURN 
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SUBROUTINE  SIMOROU ( RM , B , MOMENT ) 

C  PROGRAM  TO  PARTITION  PARTICLE  ACTIVITY  DISTRIBUTION  INTO  GROUPS 

C  USES  SIMPSON  INTEGRATION 

COMMON  RCENTR (60) ,FRAC(60) ,MC 
DIMENSION  R( 60 ) 

REAL  MOMENT , MC  ,  LRL 1 
EXTERNAL  GRANDLN 

GRANDLN  COMPUTES  NORMAL  FUNCTION  INTEGRAND  FOR  SIMPSON 
NTEGRATION 


WRITE  ( 6 , * ) *  GROUP  LEFT  RRIGHT  CENTROID  ACT 

FRAC' 

LRL  1  a  A  LOG  (  RM )  +  ( MOMENT  -  1  .  )»B**B 

C  (MUST  DECREMENT  MOMENT  TO  COMPENSATE  FOR  LOGARITHMIC  R 

INCREMENTS) 

J*0 

SUMaO . 

DP*  .  1 

DO  130  POWERa-9. , -3 • . DP 
JaJ+1 

R( J)*10. ••POWER 
RUP*10»»(P0WERfDP) 

DRaRUP-R ( J ) 

RCENTR ( J)*SQRT(R( J)»RUP) 

ARGs  R ( J ) 

F1aGRANDLN(ARG,LRL1  ,B) 

ARGa.5*(R( J)+R(J)+DR) 

F2*GRANDLN( ARG.LRL1 ,B) 

ARGaR  (  J  )  •►DR 

F3»ORANDLN( ARG.LRL1 ,B) 

DNDR  a ( 1 ./6 . )«(P1+4 .»F2  +  F3) 

FRAC ( J )  *DNDR*DR 
SUM*SUM+FRAC(J) 

C  WRITE  (6,30)  J , R ( J ) ,RUP, RCENTR(J)  ,FRAC(J) 

30  FORMAT  (13, 4(E10 . 3) ) 

130  CONTINUE 
RETURN 
END 


FUNCTION  GRANDLN(R,LRL1 ,B) 

C  COMPUTES  LOG  NORMAL  FUNC  AS  INTEORAND  OF  SIMPSON  INTEGRATION 
PROCEDURE 

REAL  LRL  1 

R2  P I =SQRT ( 2 • 3  •  14159) 

EXPARGs-. 5*(ALOG(R)-LRL1 ) ••2/B/B 
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IF (EXP ARG.LT. -200 .) EXP ARG»- 200. 
GRANDLN* 1 ./R2PI/B/R*EXP( EXPARO) 
RETURN 
END 

END  OF  FILE 
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APPENDIX  E 


Program  INTOPT  Listing 
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PROGRAM 

INTOPT1  (INPUT , OUTPUT , TAPE5* INPUT ,TAPE6 -OUTPUT ,  PLFILEaO ) 

C  COMBINES  75  EVENTS  INTO  ONE  EVENT  BASED  UPON  PARAMETER 
STATISTICS 

C  PROGRAM  COMPUTES  INTERMEDIATE  PALLOUT  FROM  ATMOSPHERIC  TEST 
SERIES  AT  NEVADA 

C  PROGRAM  USES  MODIFIED  BRIDGMAN. BIGELOW  SMEARING  EQUATION  TO 
INCLUDE  VERTICAL 

C  SPREAD  IN  CLOUD  SPACIAL  DISTRIBUTION 

COMMON  RCENTR(33) , WEIGHT( 33 ), WORK ( 1 OOOO ) 

EXTERNAL  VBL 
DIMENSION 

ACT( 150, 101 ) , ZGROUP ( 33 ) ,ATRUE(6) ,XTRUE(6) ,IPKRAT(60) 

DIMENSION  ACTCENT( 150) , ACTCUM( 1 00 , 1 00 ) , WASH ( 149) ,X( 149) 
REAL  LOONOR, NATH AN, MOMENT, MC, PI VX,P IVY 
INTEGER  EVENT 

DATA( ( ACTCUM( I , J) ,1*1 ,  100)  ,  Jal , 1 00 ) / 1 0000»0 . / 

ACTMAXsC. 

C  EVENT  DATA  POR  INTERMEDIATE  FALLOUT  COMPUTATION 

C  PARAMETERS  FOR  SINGLE  REPRESENTATIVE  BURST  FOLLOW  (YIELD 

WEIGHTED  AVGS ) 

YIELD* 1 3 • 27 
ZOKPT-34  •  1 
VECT0R*246 . 

VWIND-27.5 

C  MPH - >KM/HR 

VWIND*VWIND»1  . 60-9 
SHEAR-5. 38 

C  SET  WASHOUT  TIME  CONSTANT, HRS  (REF.  C.  E.  JUNGE,  JOURNAL  OF 
METEOROLOGY) 

WRITE(5,#) 'INPUT  TWASH  (HRS) * 

READ ( 6 , • ) TWASH 

C  SET  UP  SIZE  GROUPS,  LOG  NORMAL,  OR  LOG  NORMAL 

EXPONENTIAL  TAIL 

NGROUP*  33 

WRITE(6,») 'INPUT  NDIST  (1>L0G  NORMAL,  2>NATHANS 

W/TAIL ) ' 

READ( 5 , •) NDIST 

12  IF(NDIST.EQ.2)WRITE(6,») 'INPUT  ROLLOFF  N  AND  RM AX 

(MICRONS) • 

'  IF(NDIST.EQ.2)READ(5,*)P1 ,RMAX 
IF(NDIST.EQ.2)RMAX*RMAX»1 .E-6 
11  WRITE(6 ,•)' INPUT  RM  (MICRONS)  AND  SLOPE' 

READ(5,#)RM, SLOPE 
IF(RM.LT.O. )STOP 
RMiRM1 1 . E-6 

C  DO  1001  RM*. IE-6 , 1 .E-6 ,. IE-6 

C  DO  1000  SLOPE  =  1 .5,3-5, .25 

B» ALOG ( SLOPE ) 

M0MENT*2 . 5 

IF ( NDIST. EQ. 2) CALL  SI Z NAT ( RM , B , MOMENT , P 1 , RMAX ) 

IF ( NDIST. EQ. 1 ) CALL  SI ZCNF ( RM , B , MOMENT ) 

C  BEGIN  EVENT  LOOP  < - 

DO  710  EVENTa 1,1 


WRITE (  6  ,  •  )  *  BVENTa* .EVENT 
C  SET  OP  BURST  PARAMETERS  (MT, TOTAL  CURIES) 

TsYIELD/IOOO. 

C  MEGATONS 

PP  a  1  . 

MC*75 .  •  Y*  . 1  *  ?  F • 1 . E9 
C  MILLICURIES 

C  SET  UP  WIND  CONDITIONS  (ANGLE  IN  DEGREES ,  VELOCITY  IN 

KM/HR) 

VBCTORa-VBCTOR+270. 

C  SET  INITIAL  CLOUD  HEIGHT,  Z  DISPERSION,  YY  DISPERSION 

INXaO 
SIGY0*50. 

C  SIGYOsI . 609* 

C  *  EXP(.7+ALOG(Y)/3.-(3.25/(4.+(ALOG(Y)+5.4)**2.))) 

C  WSttG  FORMULA 

C 

ZOKFTa44  .+6. 1 «ALOG < Y) - . 205* ( ALOG ( Y >  +  2 . 42 ) • ABS ( ALOG ( Y) +2 . 42 ) 
ZOMaZOKFT/3. 28*1000. 

C  USING  HOPKINS  FORMULA 

C  LNY»ALOG(Y(EVENT)*1000.  ) 

C  Z0M*EYP(7.889+. 34*LNY+.001226*LNY**2-.005227*LNY**3 

C  *  +.000417*LNY**4) 

DO  100  K«1  ,  NGROUP 
ZGROUP ( K ) sZOM 
100  CONTINUE 

SIGZOa. 18*Z0M/10C0. 

TCa12.*ZOKFT/60.-2.5*(ZOKFT/60. )**2 
C  BEGIN  X  INCREMENT  LOOP 

DELTX-20. 

C  SET  XXsDISTANCE  DOWNWIND  PROM  GZ,  COMPUTE  ARRIVAL 

TIME 

DO  700  IN X a  1  r 149 
GRPSUMaO . 

XXaFLOAT(INX)*DELTX 

GO  TO  670 

669  PIVXaO. 

PIVYaO  . 

VECTORa38. 

XXSTARTaO. 

IP(  ..X  .  GT.  443  •  )  THEN 
VECTOR  a24 . 

XXSTARTa443 . 

PIVXa443.*COSD(38. ) 

PIVYa443.*SIND( 38. ) 

ENDIF 

IP(XX.GT. 1662. )  THEN 
VECTORa-30. 

XXST ART  a  1  6  6  2 . 

PIVXaPIVX  +  (  1662.-443. )*COSD(24  .  ) 

PIVYaPIVY+( 1662.-443.)* SIND (24.) 

ENDIF 

IF( XX .GT . 2069 . )  THEN 


VECTORa-MO . 

XXSTART*2069 • 

PIVXaPIVX+(2069.-l662. )»C0SD(-30.) 

PI VYaPIVY+( 20  69. -1662  . ) *SIND ( -  30 . ) 

ENDIF 

670  TAaXX/VWIND 

DELTATaDELTX/VWIND 
DELSBC»DELTAT*3600 . 

C  COMPOTE  CLOOD  DISPERSION  IN  YY  DIRECTION  (KM) 

TSTARaMIN ( 3 • ,TA) 

C  SIQYaSQRT 

C 

( SIOYO**2 . • ( 1 .+(8.»TSTAR/TC))+.5#(SHEAR»SIGZ0»TA)»»2. ) 

SIOIa50.+555.*(EXP(-TA/MO. )-EXP (-TA/ 10 . ) ) 

C  WRITE ( 6  ,  • )  '  ZOMa ' , ZOM, '  SIGYa • , SIGY , *  TAa'.TA 

C  COMPOTE  ALTITODB  OP  EACH  CONTRIBOTING  PARTICLE 

SIZE  GROOP 

DO  680  Kal  ,  NGROUP 
RPaRCBNTR ( K ) 

IF(ZGROOP(K).GT.O.)THEN 

ZBaZGROOP(K) 

ELSE  ZBaO. 

ENDIP 

DZaVEL(ZBfRP)*DBLSBC 

IP( (ZGROOP(K) .GT.O.  ) . AND . 

((ZGROUP(K)-DZ) .LT.O.  )) 

*  ‘  RDEPaRCENTR(K) 

ZGROOP(K)aZGROOP(K)-DZ 

GRPSUMaGRPSUM* 

GRNDF(ZOROOP(K) ,SIGZO,TAfO. )• WEIGHT (K) 

*  *MC 

C  WRITE(6,679)INX,K,ZB, DZ, 

C  *  GRNDP(ZGROUP(K) ,SIGZ0,TA,0. ) ,GRPSUM 

679  FORMAT (  •  INX',13,'  GP',13,'  ZB',E10.3,' 

DZ' , B 1 0 . 3 i 

*  '  GRNDP ' , E 1 C . 3  » '  GRPSOM*  ,E10.3) 

680  CONTINOE 

C  SET  VALOE  OF  CUMULATIVE  ACTIVITY  ON 

CENTER LINE (MILL ICURIES) 

SEDaGRPSUM* EXP (-TA/T WASH) 

WASH ( INX )aMC»( 1 . -EXP ( -TA/TWASH ) ) 

ACTCENT( INX) aWASH( INX)+SED 

C  WRITE(6,685)XX,TA,SED,WASH(INX) , RDEP 

685  FORMAT (  '  XXa',E8.3,'  TAa',E8.3,'  SEDa',E8.3 

*  ,»  WASHa' ,E8. 3, '  RDEPa'fE8.3) 

C  WRITE(6, •) ' INXs' ,INX, '  ACTCENT a ' , ACTCENT ( INX ) 

GO  TO  390 

C  COMPUTE  ACTIVITY/AREA  ON  THE  GROUND  ('ACT'  IN 

CURIES/SQ  M) 

C  FOR  EACH  XX  POSITION,  VARY  YY  OUT  TO  +-100  KM 

DO  380  Jal ,  101 
YYaFLOAT( J-1)»20.-1000. 

EXP0Na.5»( YY/SIGY)»»2. 

IF  (EXPON.GT.  100.  )  THEN 
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BXPFACaO.O 

ELSE 

SXPFACaEXP(-EXPON) 

ENDXF 

FYTA« (SQRT( 2. •3.141 59) #SIGY )••(-! . )»EXPFAC 
IF(INX.EQ.I)  THEN 

ACT(INX, J)aACTCENT(INX) 

*  •FYTA/DELTX*2.59 

C  MILLICURIBS  PER  SQOARE  MILE 

ELSE 

'  ‘  ACT(INX, J)a(ACTCENT(INX) 

A  -ACTCENT( INX-1 ) )»FYTA/DELTX 

BNDIF 

C  UPDATE  CUMULATIVE  ACTIVITY 

XCapIVX+(XX-XXSTART)*COSD( VECTOR) 

*  -YY*SIND( VECTOR) 

YC*PIVY+(XX-XXSTART) •SIND (VECTOR) 

A  +YI»COSD( VECTOR) 

ICaINT( (XC  +  30  00.  )/60. )-2000./60. 
IF(IC.LE.O)  QO  TO  690 
JCaINT((YC+3000.)/60.) 

ACTCUM(IC,JC)aACTCUM(IC, JC)+ACT(INX, J)/9. 


,  ACT ( INX , J ) 


ACTMAXaMAX (ACTMAX,ACTCUM(IC,JC) ) 

690  CONTINUE 

C  WRITE ( 6 , • ) *  XXa*,XX»'  YYa  *  ,  YY , '  ACTa' 

380  CONTINUE 

390  CONTINUE 

700  CONTINUE 

710  CONTINUE 

WRITE(6,#) ' ACTMAX  = ' , ACTMAX 

GO  TO  212 

C  COMPUTE  ACTIVITY  ENCLOSED  BY  EACH  CONTOUR 
C  EACH  60  X  60  KM  CELL  IS  1390  SQUARE  MILES 
208  DATA  DOWNIO .DOWN20 , D0WN30 .DOWN40 ,D0WN50 .D0WN60  /6»0./ 

AREA  1 0aO . 

DO  200  I a  1 , ICO 

DO  210  J= 1 , 100 

IF ( ACTCUM (I ,  J ) .GT. 1 0 . ) AREA  1 0  a AR E A  1 0  + 1 390 . 

IF(ACTCUM(I,J) .GT.60.)  THEN 

DOWN60aD0WN60  + 1390  .  •ACTCUMd, J) 

ACTCUMd,  J)a60. 

ENDIF 

IF( (ACTCUM(I, J) .GT. 50 . ) .AND. ( ACT  CUM ( I , J ) . LE . 60 . ) ) 

A  DOWN  50  a  DOWNS  0+1  390  .  •ACTCUMd  ,  J) 

IP  (  (ACTCUMd,  J)  .GT.  40.  )  .AND.  (ACTCUMd,  J)  .LE.50.  )  ) 

*  DOWN  40  =  DOWN40  +  1  390  .•ACTCUMd,  J) 

IF  (  (  ACTCUMd,  J)  .GT  .  30  .  )  .AND.  (  ACTCUM  ( I ,  J  )  .  LE  .  40  .  )  ) 

*  DOWN 30 a DOWN 30+ 1  390. • ACTCUM ( I , J ) 

IF  (  (  ACTCUMd,  J)  .GT.  20  .  )  .  AND.  (  ACTCUM  (  I ,  J  )  .  LE  .  30  .  )  ) 

*D0WN20aD0WN20+1 390  .  •  ACTCUMd  ,  J) 

I F  (  (  ACTCUMd  ,  J)  .GT  .  10  .  )  .AND.  (  ACTCUMd  ,  J)  .LE.20  .  )  ) 

*  DOWN  lOaDOWN  10+1  ?90  .  •ACTCUMd  ,  J) 

210  CONTINUE 
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200  CONTINUE 

CUM10*DOWN60+UOWN50*DOWN4O4-DOWN  30+D0WN20+DOWN 1 0 
WRITE ( 6 , 209) D0WH60 , D0WN50 , DOWN40 , DOWN  30 ,DOWN20 , DOWN  1 0 
209  FORMAK  •  D0WN60-  1 0  =  •  ,  6(  E 1  0 . 3 )  ) 

WRITE ( 6  ,  •  )  'CUMULATIVE  DOWN  UP  TO  CONTOUR  10=',C0M10 
WRITE(6 ,•) 'CUMULATIVE  AREA  UP  TO  CONTOUR  10=',AREA10 
212  CONTINUE 

C  COMPUTE  FIGURE  OF  MERIT  FOR  SIZE  DISTRIBUTION  OPTIMIZATION 
DATA 

( ATRUEC I ) i I s 1 ,6)/0. , 14.7E6,47.1E6,55.2E6,60.3E6,63.3E6/ 

DATA(XTRUE(I) , I* 1 , 6 ) /C . t 480 . , 1 700 . , 2 1 00 . , 2460 . , 2800 . / 
F0  =  0  . 

FOa( ATRUE( 2) -ACTCENT (24) )«*2  . 
F0sF0*(ATRUE(3)*ACTCENT(85))**2. 

F0«P0+(ATRUE( 4) -ACTCENT ( 105) )»»2. 
P0=F0+(ATRUB(5)-ACTCENT( 123) )#t2. 
F0=F0+(ATRUE(6)-ACTCENT( 140) )»»2. 

WRITE(6,#)ATRUE(2) , ACTCENT(24) 

WRITE(6,*)ATRUE(3) .ACTCENT (85) 

WRITE(6,«>)ATRUE(4)  ,  ACTCENT(  105) 

WRITE(6 ,*)ATRUE(5) , ACTCENT (123) 

WRITE(6,*)ATRUE(6) , ACTCENT ( 1 40 ) 

WRITE ( 6 , • ) ' RM  =  ' , RM , '  SLOPE* ', SLOPE  ,  ’  F0=’,F0 

1000  CONTINUE 

1001  GO  TO  11 
GO  TO  800 

C  DISSPLA  CONTOUR  PLOTTING  CALLS  FOLLOW - 

C  PACKAGE  PLOTS  ACTIVITY  GROUNDED  VS  TIME 
799  CALL  COMPRS 

CALL  P AGE ( 8 . 5  i  1  1  •  ) 

CALL  AREA2D(5. ,5.  ) 

CALL  LINES('DATAffIpoAYt1) 

CALL  LINES( 'WASHOUT  ONLY,  IpKRAy  2) 

CALL  LINES (  *  WITH  SEDIMENT,’  IpRRAY  3) 

CALL  HEADIN( ' ACTIVITY  GROUNDED  VS  DISTANCE,  10Q  1  ?5  1 
CALL  INTAXS 

CALL  XNAME( 'KILOMETERS,  1Q0) 

CALL  YNAME( 'KILOCURIES, ’ 1Q0) 

CALL  OR AF (0. ,5 00. ,3000. ,0., 10. ,100.) 

ENCODE ( 100 , 15 .LABEL )NDIST ,RM, SLOPE 
15  FORMAT ( '  DIST',12,',  RM=',E8.3>',  SLOPEs ',  E8 . 3  } 

CALL  RLMESS(L ABEL,  1 00 , 1  00  .  ,75  .  ) 

CALL  THKFRM ( . 02 ) 

CALL  FRAME 
DO  10  1=1,149 
X(I) =20.*FL0AT(I) 

WASH(I)=WASH(I)/1 . E6 
ACTCENT( I ) =ACTCENT( I ) / 1 .E6 
C  MILLICURIES-->KILOCURIES 

10  CONTINUE 
DO  20  1=1,6 


ATRUE(I)*ATRUE(I)/1  .  E6 
20  CONTINUE 

CALL  MARKER(  15) 

CALL  CURVE(XTRUE,ATRUE,6,-1  ) 
CALL  MARKER( 5) 

CALL  CURVE(X,WASH, 149 , 15) 
CALL  MARKER (  17) 

CALL  CURVB(X,ACTCENT, 149 , 15) 
CALL  LEGBND(IPKRAY,3, .25,4.  ) 
CALL  BNDPL(O) 

CALL  DONEPL 
STOP 

800  CONTINUE 
END 

CKMIMMIMMfllHMMtlMIMttl 
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SUBROUTINE  SIZCNF ( RM , B , MOMENT ) 

C  SUBROUTINE  TO  PARTITION  PARTICLE  MASS  DISTRIBUTIONS  INTO 

GROUPS 

COMMON  RCENTR(33) ,MASPRA(33) 

DIMENSION  RLEFT(33) ,RRIGHT(33) ,MASCNF(33) 

REAL  MASCNF,MASFRA,MOMENT,MC 
INTEGER  GROUP, OFFSET 
EXTERNAL  CNF 

C  CNF  IS  THE  CUMULATIVE  NORMAL 

FUNCTION ( RADIUS ,RM, BETA, MOMENT) 

C  DEFINE  LOG  NORMAL  FUNCTION  PARAMETERS,  RM  BETA 

3C  FORMAT  (I3,5(E10.3)  ) 

C  WRITE(6,»)  *  RM=  ' , RM ,  •  B=',B,*  MOMENT =', MOMENT 

WRITE(6,») • PURE  CUMULATIVE  NORMAL  SIZE  DISTRIBUTION 
OPERATIVE* 

C  WRITE (6,40) 

40  FORMAT ( *  GROUP  RLEFT  RCENTROID  RR IGHTOID  ACTV  FRAC 
CUMULA* ) 

EXPONs-6 . 

GROUP* 1 

RLEFT(GR0UP)*0. 

RRIGHT( GROUP )= 10 .••EXPON 
RCENTR ( GROUP ) *RRIGHT( GROUP ) 

C  WRITE ( 6 , 60 )  RRIGHT(GROUP) , RM, B 
C60  FORMAT  (  3E10.3) 

MASCNF( GROUP) = CNF (RRIGHTC GROUP) , RM,B, MOMENT) 

MASFR A ( GROUP )*MASCNF( GROUP) 

DO  130  GROUP  *  2 ,  32 
EXPON*  EXPON-*- .  1 
RRIGHT( GROUP) =10 EXPON 
RLEFT (GROUP) *RRIGHT( GROUP- 1 ) 

RCENTR ( GROUP )=SQRT( RLEFT ( GROUP )• BRIGHT ( GROUP ) ) 
R*RRIGHT( GROUP) 

MASCNF( GROUP) =CNF(R , RM,B , MOMENT) 
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MASFR A ( GROUP ) eMASCNF ( GROUP ) -MASCNF ( GROUP-  1  ) 

130  CONTINUE 

C  LUMP  LARGE  SIZES  TOGETHER  IN  GROUP  WITH  RCENTRslOOO  MICRONS 

K*  33 

RRIGHT(K) «999. 

RLEFT(K)»RRIGHT(K-1) 

RCENTR(K)-RRIOHT( K-1  ) 

MASFRA(K)a1 . -MASCNF ( K  -  1  ) 

MASCNF(K)*CNF(RRIGHT(K) , RM ,B .MOMENT ) 

C  OUTPUT  GROUP  DIVISIONS 
DO  140  X* 1 , 33 
C 

WRITE(6,30)K,RLEFT(K) , RCENTR ( K ) , RR IGH * v  a. ) , MASFR A ( K ) .MASCNF ( K ) 

140  CONTINUE 
RETURN 
END 

CIIMIMIMtMIIIIIIIMIIMIIIIIIMIIIttlltllllllttKIIMIItlKtll 
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SUBROUTINE  SI ZN AT ( RM , B , MOMENT , P 1 , RMAX ) 

COMMON  RCENTR( 33) , ACTFRA( 33) 

C  SUBROUTINE  PARTITIONS  NATHANS  R4  DISTRIBUTION  INTO  ACTIVITY 

GROUPS 

DIMENSION  RLEFT(33) ,RRIGHT(33) , ACTCUM(33) 

REAL  MOMENT5MC,IA,IB 
INTEGER  GROUP, OPFSET 
EXTERNAL  CNF 

CNF  IS  THE  CUMU'/.IVE  NORMAL  PUNCTION 
DEFINE  NATHANS  PARAMETERS 

RO»RM«EXP((P1-1  .  ) • B*B ) 

R2PIsSQRT( 2 . *3 .  14159) 

30  FORMAT  (I3»5(E10.3)) 

C  WRITE ( 6,40) 

40  FORMAT ( '  GROUP  RLEFT  RRIGHT  RCENTR  ACTV  FRAC 

CUM'  ) 

EXPONn-6 .  1 
G  ROUPs  0 

60  FORMAT  (  3E10.3) 

DO  280  GROUPsI  ,  32 
EXPON »  EXPON*  .  1 
RRIGHT(GROUP)»10.##EXPON 
IF  (GROUP. EQ.1)  THEN 
RLEFT(GROUP) aO. 

RCENTR ( GROUP )* RR IGHT ( GROUP  ) 

ELSE 

RLEFT (GROUP) a RRIGHT (CROUP-  1  ) 

RCENTR(  GROUP  )»SORT(  RLFFT(G1  0'JP)»RR  IGHT(  GROUP  )  ) 

END  IF 

RP=R  RIGHT  (0*.  ■  ; 

C  CALCULATE  D  F  N  0  M  ±  N  «*  'JR  WHICH  I?  P  0 i C p ‘  I  yK  4L  1 :  TOTAL 

ACTIVITY 

IA*R0**(P1/2.4’.5),RM##(  5-M  MENT)*B*R2P 
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EXP ( .5*M0MENT»M0MENT»B»B)*CNF( RO , RM,B, MOMENT) 

IF  (MOMENT. EQ. (PI-1.))  THEN 
IBaRO**P1*ALOG(RMAX/RO) 

ELSE 

IEaRO**P1*(RMAX**(  MOMENT-P  1  1  .  ) -RO*»  ( MOMENT- P  1  ♦  1  . 

)) 

*  /(MOMENT-PI + 1  .  ) 

END  IF 

DBNOMalA  4-  IB 
FaIA/(IA  ♦  IB) 

C  NOW  CALCULATE  NUMERATOR  WHICH  IS  PROPORTIONAL  TO  ACTIVITY 

DOWN 

IF  ( RP . LT  .  RO ) THEN 

r»ISO»BT.(,Ue^86  DU  O  UnUDUT^ 

\  A*  *  |  till  |  W  |  UViiWkl  *  / 

ELSE 

CNPPARTsCNP(RO ,RM,B, MOMENT) 

END  IF 

VI A»R0** ( .5+P1/2 . )»RM»»( . 5-P1/2. ) »RM«» ( MOMENT ) »B» R2 

PI» 

EXP (MOMENT *M0MENT*B*B* ,5)*CNFPART 
IP  ( RP . LT  .  RO ) THEN 
VIB=0 . 

00  TO  250 
END  IF 

IF  (MOMENT. EQ. (PI-1 . ))THEN 
VIB=RO*»P1»ALOG(RP/RO) 

ELSE 

VIBaRO*#Pl*(RP#,(  MOMENT-P  1  ♦  1  .) -R0»»  (  MOMENT-P  1  •>  1  .  ) 

) 

/ ( MOMENT-P 1  ♦  1  .  ) 

END  IF 

250  CONTINUE 

ACTCUM ( GROUP) s ( V I A  +  V IB ) / DENOM 
IF  (GTOUP. LQ. 1 )  THEN 

ACTFRA ( GROUP )« ACTCUM (GROUP) 

ELSE 

ACTFRA ( GROUP )= ACTCUM (GROUP )- ACTCUM (GROUP- 1 ) 

END  IF 

C  WRITE (6,30) 

GROUP, RLEFT( GROUP ) ,RRIOHT( GROUP ) , RCENTR( GROUP) 

C  *  ,  ACTFRA ( GROUP ) , ACTCUM ( G ROUP ) 

280  CONTINUE 

C  DEFINE  GROUP  WHICH  LUMPS  BIG  PARTICLES  TOGETHER  WITH 

RCENTRa 1000  MICRONS 
GRCUPs  33 

F.L EFT  (  0 ROUP  )*RRIGHT(GROU P-1  ) 

RRIOHT(OROUP) sRMAX 
RCENTR(GR0UP)sRRIGHT(GRCUF-1) 

ACTFRA( GROUP) a  1 . -ACTCUM' GROUP- 1 ) 

ACTCUM (GROUP) a ACTCUM ( GROUP- 1 )+ACTFR A ( GROUP) 

WRITE ( 6 , 30 ) 

GRO UP, RLEFT( GROUP ) , RRIGHT ( GROU? ) , RCENTR (GROUP) 

“  ,  ACTFRA(GROUP) , ACTCUM ( GROUP ) 
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RETURN 

END 

CIMMIMIM 


c 


FUNCTION  CNF(R ,RM,B, MOMENT) 

C  COMPUTES  CUMULATIVE  NORMAL  FUNCTION  FOR  ACTIVITY 

DISTRIBUTION 

REAL  MOMENT 
R3»RM*EXP(MOMENT*B*B) 

Z*( ALOG(R)-ALOO(R3) )/B 
C  WRITE( 6,50)  Z 

C  50  FORMAT  ( *  Z»' , E 1 0  -  3  > 

IF  (Z)  1060,1050,1050 

1050  CNF*1 .-.5/( 1 196854»Z*. 1 1 5 1 94*Z««2  ♦ . 000344*Z»*3 

♦.019527*Z 

*  •  •  14  )  #  •  4 


GO  TO  1070 
1060  Za-Z 

CNF  = .5/(1 196854*Z+. 115194* Z ••2^.000 344* Z»»3+. 01 95 27 

•  Z  *  •  4  ) 

~  ««4 


1070  RETURN 
END 


C 


C 


FUNCTION  GRNDF(Z3,SIGZ0,TA,ZTR0P) 

REAL  KDZ 
KDZaO. 

SIGZsSIGZO* 1000.>SQRT(2.*KDZ*TA*3600. ) 

DELTZsZB-ZTROP 

AsOELTZ/SIGZ 

IF  (A)  1060,1050,1050 

1050  CUMNFsl .-.5/( 1 .+.196854»A+. 115194*A»»2  + . 000344*A«*3 
.  0 1  9527*  A 

*  *«4  )«»4 

GO  TO  1070 
1060  As-A 

CUMNFa  . 5/(1. ♦.196854»A+.115194*A#»2-*..0C'0344»A»B3+. 01952 

7»A«*4) 

*  »§4 

1070  GRNDFa 1 . -CUMNF 
RETURN 
END 


FUNCTION  VFL ( ZB , RP ) 

C  COMPUTES  TERMINAL  VELOCITY  OF  SPHERICAL  PARTICLES  USING 

FULL  DELFIC 
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nwjiuw*m’ 


C  FALL  MECHANICS  AND  O.S.  STANDARD  ATMOSPHERE  EQUATIONS 

C  SUBROUTINE  ATMOS  RETURNS  DENSITY,  VISCOSITY  AND  MEAN  FREE 

PATH  AT  ALT  ZB 

REAL  MPP , LOOREY 
C  ' EXTERNAL  ATMOS 

CALL  ATMOS(ZB, DENS, VISC, MFP) 

DENSF*2600. 

0*9.8 

QZa4*DENS*(DBNSP-DENS)*G*(2.*RP)**3/(3.*VISC*VISC) 

C  WRITE(6,1051)QZ, DENS, VISC, MFP 

Cl 051  FORMAT ( * QZ, DENS, VISC, MFP » , 4(E10 .3)  ) 

YY* ALOQ ( QZ ) 

C  DENSF  IS  PARTICLE  DENSITY,  Q  IS  ORAV ,  QZ  IS  DAVIES  NUMBER 

DEFINED  IN 

C  DNA  TR-5159F-1 ,  P24 

IF  ( QZ  .LE.  .326DTHEN 

VEL*VISC«QZ/(24 . »DENS*2 .  *RP) 

00  TO  1100 
ENDIF 

IF ( ( QZ  .GT.  .3261)  .AND.  (QZ  .LE.  84.175))THEN 

VEL*(VISC/(DENS»2. *RP ) )*EXP(-3. 1 8657+ . 992696* YY- 

* 

1 .53193E-3*YY*»2-9.87059E-4*YY**3-5.78878E-4*YY**4 

*  ♦8.551759E-5*YY»*5-3.27815E-6*YY**6) 

GO  TO  1100 

ENDIF 

IF ( ( QZ  .GT.  84.175)  .AND.  ( QZ .LT . 1 40 ) ) THEN 

VEL*VISC*QZ*(4. 1 66667E-2-2 . 3363E-4*QZ+2 . 0 1 54E-6*Q 

Z*  *2 

-6.9105E-9*QZ**3)/(DENS*2.*RP) 

GO  TO  1100 
ENDIF 

IF ( ( QZ  . GE .  140)  .AND.  ( QZ  .LT.  4.5E7))THEN 

LOGREYa-1 . 295 36^. 986* YY/2. 303-. 046677* (YY/2 . 303) * 

•2 

*  +.001 1235*(YY/2.303)**3 
REYN0L*10.**L0GREY 
VEL*REYN0L*VISC/(2. *DENS*RP ) 

ENDIF 

C  KNUDSEN  SLIP  CORRECTION 

1100  IP(RP.CT. 100.»MFP)  THEN 

EXPFACaO. 

ELSE 

EXPF AC* EXP (-. 656*2 . *RP/MFP ) 

ENDIF 

VEL  *VEL* ( 1 .  ♦  (  1 ,644  +  .552*EXPFAC)*MFP/ 

*  ( 2 . • RP  )  ) 

RETURN 

END 


C****»*«» 


SUBROUTINE  ATMOS ( ZB , DENS, VISC, MFP) 


C  SUBROUTINE  CALCULATES  ATMOSPHERIC  TEMPERATURE ,  PRESSURE  AND 

•  -VISCOSITY 

C  FOR  ALTITUDES  UP  TO  84,852  METERS  (MRS  UNITS  EMPLOYED) 

REAL  LX , NUMDEN , MPP 
I  .  ^ IF (ZB  .LB.  11000.) THEN 

^LK*-. 006545 
TK»288 . 15 
PK*  101 300 . 

PB*PK* ( TK/TB ) • • ( .034164/LK) 

....  •  00  TO  1240 

ENDIF 

IF ( ( ZB  .OT.  11000.)  .AND.  (ZB  .LE.  20000. ))THEN 
ZKallOOO. 

mtr  -  4  C  c  p 

.  i  tv*c  1  W  • 

?X=2269G. 

LKsO. 

TB*TK 

PB*PK»EXP(-.034l64«(ZB-ZK)/TK) 

00  TO  1240 
ENDIF 

IF ( (ZB  .GT.  20000.)  .AND.  (ZB  .LE.  32000. ))THEN 
ZK*20000. 

TK-216.65 

PK.5528. 

LK« . 00 1 

TB«TK+LK»(ZB-ZK) 

PBaPK»(TK/TB)«*( .034164 /LK ) 

00  TO  1240 
ENDIF 

IF ( ( ZB  , GT .  32000.)  .AND.  (ZB  .LE.  47000. ))THEN 
ZK*  32  000 . 

TKs  228 . 65 
PK*888 . 8 
LKs.0028 

TBoTK+LK»( Z3-ZK) 

PB*PK»(TK/TB) »•( .034164/LK) 

GO  TO  1240 
ENDIF 

IF ( ( ZB  .GT.  47000.)  .AND.  (ZB  .LE.  51000. ))THEN 
ZK.47000. 

TKr 270 . 65 
PKa 1 10.873 
LK*0  . 

TB*TK 

PB*PK»EXP(-.034164»(ZB-ZK)/TK) 

GO  TO  1240 
ENDIF 

IF ( ( Z3  .GT.  51000.)  .AND.  (ZB  .LE.  71000. ))THEN 
ZK»5 1 000 . 

T  K  *  2  7  0 . 65 
PKs66 .9218 
LK  =  -  . 0028 
TB*TK+LKf(ZB-ZK) 


or  n 


PB»PK*(TK/TB)«*( .034164/LK) 

00  TO  1240 
ENDIF 

IF ( ( ZB  ,0T.  71000.)  .AND.  (ZB  .LE.  84852. ))THEN 
ZK»7 1 000 . 

TK.214 . 65 
PK*3. 95536 
LX»- . 002 

TB*TK+LK# ( ZB-ZK  ) 

PB*PK»(TK/TB)»»( .034164/LK) 

00  TO  1240 
ENDIF 

IF ( ( ZB  .OT.  84852 ) ) THEN 
WRITE( 6 ,1230) 

00  TO  1241 
ENDIF 

1230  F0RMAT( 'PARTICLE  ALTITUDE  BEYOND  RANGE  OF  PROGRAM 
APPLICABILITY' ) 

1240  VISCs 1 .46E-6«TB»»1 . 5 /  (  TB+ 1 10.4) 

DENS«.003484»PB/TB 
NUMDEN«DENS»2.55E25/1 .23 

PER  CUBIC  METER:  2.55E25  MOLECULES  PER  CUBE  AT  SEA 

EVEL. 

DENSITY  AT  SEA  LEVEL  IS  1.23  KG/CUBE 
MFP  *  1 ./(I .414»NUMDEN»4.5E-19) 

C  MFP*  1 ./(R00T2*N*SIGMA) 

1241  RETURN 
END 

END  OP  FILE 
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